ATR trajectory tracking system (A-Track)

ABSTRACT

A-Track is a tracking process that is driven by automatic target recognition techniques (ATR), mobility analysis using digital maps and exploitation of other constraining information. It is structured as an optimization problem amenable to the formalism from classical relaxation labeling algorithm. The novel combination of techniques in A-Track provides a new approach to the temporal problem of establishing and extending target tracks through a sliding time window involving a sequence of multiple data frames. A data frame is a collection of sensor reports taken from one scan of a predetermined surveillance area that is being systematically and repeatedly scanned over time. The relative time differences between reports within a frame are generally much smaller than the time difference between different frames. A-Track is especially designed to handle cases having relatively large time intervals between frames that cause problems for conventional tracking algorithms based on predictive models.

CROSS REFERENCE TO RELATED APPLICATION

This application is a continuation-in-part of U.S. patent applicationSer. No. 10/370,307, filed 10 Feb. 2003 now abandoned, which is herebyincorporated by reference.

RIGHTS OF THE GOVERNMENT

The invention described herein may be manufactured and used by or forthe government of the United States for all governmental purposeswithout the payment of any royalty.

FIELD OF THE INVENTION

The invention relates generally to computerized techniques forprocessing data obtained from target recognition sensors obtained totrack multiple discrete objects.

BACKGROUND OF THE INVENTION

Remote surveillance of airborne and ground vehicles rely upon strategicor theatre-level systems (e.g., satellites, AWACS) for detecting andidentifying targets by sensing radiated or reflected electromagneticenergy from target vehicles. Generally, these surveillance systemsprovide automatic target recognition (ATR) capabilities wherein themulti-spectral and spatial characteristics of each return can identifythe type of vehicle. Enhanced range or surveillance or additionalspectral bandwidth is often obtained by data fusion, using data from anumber of different surveillance platforms.

Once one or more targets of interest have been located and identified,cooperative action may be directed against the targets. The location andcharacterization of the target of interest is handed off to a tacticalplatform (e.g., aircraft, missile, directed energy weapon) for action.These tactical platforms require a high-degree of knowledge of currentand predicted future location of the target in order to engage. Thisresolution of tracking is not available from the surveillance platforms.Consequently, a tracking system (e.g., tracking radar system) rapidlyscans a narrow aerial volume or surface area about the location receivedfrom the surveillance system. The repetition rate of sweeps producesframes of returns that have time difference between respective framesthat is relatively short as compared to the mobility of the identifiedtargets. Thus, linking returns between frames is a straightforwardmatter of linking the most closely spatially related returns betweenframes. Predicting future locations of the target is made byextrapolating the track of previous returns.

Even with rapid sweep repetitions, the amount of processing required ishigh. Attempts to increase the number of targets that may be trackedgenerally requires increasing the time difference between returns, andthus adding a degree of complexity to the predictions required to linkreturns. For instance, U.S. Pat. No. 5,537,119 to Poore, Jr. is anexample of a predictive technique in tracking targets. In particular,Poore, Jr. relies upon an iterative Lagrangian Relaxation technique tolink returns between each processed frame in order to assign a trackrespectively to each target.

While Poore, Jr. provides certain processing advantages over other knownpreditive techniques, Poore, Jr. shares the general limitation of otherpredictive techniques in that time differences between frames must stillbe relatively short as compared to the mobility of the targets. Inaddition, these predictive techniques presume that each target ispresent frame to frame in order to maintain a track. Often, returns fora target are not present in a given frame due to factors such as aspectsof the target that present a small return, relative velocities,intervening obstructions or atmospheric interference, limitations of thetargeting system, etc.

Consequently, a significant need exists for a system and method fortracking that allows larger time intervals between frames and that maymore robustly handle tracking when returns are missing in certainframes.

SUMMARY OF THE INVENTION

The present invention addresses these and other problems in the priorart with a method and system that can use Automatic Target Recognition(ATR) techniques, even with large temporal spacing between targetreports, to develop a trajectory description output. Moreover, a robustapproach is described that accommodates complications associated withlarge temporal spacing.

In one aspect of the invention, the system and method determines thetrajectories of a plurality of periodically sensed targets byiteratively grouping received sensor reports into a frame, linking eachreport in the frame to a previous frame with an initial assignmentvector, converting the linking matrix into a plurality of compatibilitymatrices, and optimizing an objective function formed from thecompatibility matrices such that the assignment vectors converge to atrajectory description output. Thereby, large temporal spacing isaccommodated, even with a significant likelihood that any given reportmay be associated with a plurality of trajectory tracks being monitored.

In another aspect of the invention, optimizing an objective functionthat is advantageously utilizes a relaxation label algorithm (RLA) loop.In particular, a support function is defined for each compatibilitymatrix as a weighted average of the respective compatibility matrix withthe assignment vector. The weighted sum of the support functions isassembled into an objective function. Then the RLA loop is iterativelyperformed to optimize the gradient of the objective function, using theupdated assignment vector each iteration. Ideally, the RLA loop willconverge to a trajectory description output for the plurality ofperiodically sensed targets.

The marriage of ATR and tracking techniques provides robust tracking insituations where a dedicated tracking system is not available. Forexample, situations where large volumes of terrain or airspace are beingmonitored, slow scan rates are being used, and/or many closely spacedtargets are being monitored create situations where reports areambiguous. An “A-track system” is presented that works through theprobable assignments of for each report to determine the trajectorytracks represented by the reports.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a diagram of target reports received from a plurality ofsurveillance windows within a single time frame.

FIG. 2 is a depiction of a plurality of report frames in a sliding timewindow, with reports forming tracks across subsequent frames.

FIG. 3 is a diagram of a missing label determination based on mobilityanalysis.

FIG. 4 is a diagram of determining a probability of an observed reportbased on mobility analysis.

FIG. 5 is a diagram of a tracking procedure for tracking the reports ofFIG. 2.

FIGS. 6–7 together form a more detailed implementation of the trackingprocedure of FIG. 5.

The accompanying drawings, which are incorporated in and constitute apart of this specification, illustrate embodiments of the invention,and, together with the general description of the invention given above,and the detailed description of the embodiments given below, serve toexplain the principles of the present invention.

DETAILED DESCRIPTION OF THE INVENTION

1. Introduction

A-Track is a tracking process that is driven by automatic targetrecognition techniques (ATR) and structured as an optimization problemamenable to the classical relaxation labeling algorithm. The novelcombination of techniques in A-Track provides a new approach to thetemporal problem of establishing and extending target tracks through asliding time window involving a sequence of multiple data frames. A dataframe is a collection of sensor reports taken from one scan of apredetermined surveillance area that is being systematically andrepeatedly scanned over time. The relative time differences betweenreports within a frame are generally much smaller than the timedifference between different frames. Tracking based on high rangeresolution radar (HRR) sensor reports from ground targets is theapplication used to describe the A-Track approach; however, the approachis applicable and extendable for tracking problems based on othersensors.

The ‘A’ in A-Track emphasizes that it is an ATR-driven tracking methodthat integrates the information from sensor signatures into the problemof associating reports across different frames. This is especiallysignificant when the time intervals between frames are relatively large;i.e., large enough to cause problems for conventional tracking methodsbased on predictive models. Over large time intervals, there are nosimple choices for selecting dynamic models to predict ground targettrajectories. Basic physical constraints and behavioral heuristics areused to establish a set of special matrices that are functions of thegiven sensor data and available knowledge bases.

With reference to FIG. 1, an illustrative depiction of an A-Track system10 receiving a frame input 12 from a plurality of sensors 14–18, eachproviding a portion of a plurality of reports 20 from a respectivesurveillance area 24–28. The A-Track system 10 processes the frameinputs 12 to generate a trajectory report 30.

Analysis used in the A-Track system 10 is based on establishing therelative likelihood of reports (α,β,χ,δ) in different frames (Ω),depicted in FIG. 2, being associated with the same target or track index(a, b, c, i). The track solution is specified by a set of distributionsthat stipulates the relative weights for labeling track segments withina given frame with the sensor report labels associated with that frame.This set of report labels also includes a null label that is associatedwith missing reports in a frame for some of the target tracks. A-Trackemploys an enhanced version of the relaxation algorithm tosystematically organize and integrate available information into anoptimization problem that can be solved with a simplified iterativeupdate rule. The iterative relaxation process takes place in a slidingtime window of F+1 frames that contains a fixed number of most recentlycollected frames. The relaxation algorithm modifies only variables inthe active window.

The reports and tracks depicted in FIG. 2 illustrate several problemsovercome by the A-track system 10. For example, the time spacing may belarge enough between frames that kinematic, predictive techniques areinadequate. The time spacing makes several reports a possible candidatefor assigning to a track, as illustrated by the dotted line portions ofthe tracks that are possible tracks along with the solid line of theactual track that has to be determined.

The reports and tracks depicted in FIG. 2 further illustratecomplicating factors that are addressed by the A-track system 10.Occasionally, a report is not received on a target for which a track hasbeen assigned, or a report is sensed for a target for which no track hasbeen previously assigned. The report may be missed due to ATRlimitations, such as interference, target aspect angle, masking terrain,or a target that has no relative movement with respect to the sensor orterrain.

Below, a brief overview of standard formalism and conventional notationgenerally known for relaxation labeling algorithms (RLA) is describedunder subheading “2. Relaxation Labeling Algorithm”. Then, theconventional notation is customized to make it suitable for A-Trackunder the subheading “3. Modified RLA for A-Track”. Here severaloriginal theorems are presented to enhance the basic RLA and methods areoutlined to deal with a multiplicity of compatibility matrices thatfacilitate the integration of diversified modeling constraints. The mostbasic and complex compatibility matrix, which is utilized in multipleways throughout the A-Track model, is developed in Section 4. Sections 5and 8 present two additional compatibility matrices that complete theconstraint modeling which encompasses the contextual information for theHRR-driven tracking problem. Section 6 addresses the methodology forintegrating multiple support functions and Section 7 outlines algorithmsfor initiating new tracks and removing lost tracks. Section 9 provides asummary of the features incorporated in A-Track.

2. Relaxation Labeling Algorithm

Relaxation labeling algorithms (RLA) are non-linear and paralleliterative procedures introduced by Rosenfeld, Hummel, and Zucker.(Rosenfeld A., R. A. Hummel, and S. W. Zucker, ‘Sceen labeling byrelaxation operations,’ IEE Trans. Syst. Man Cyber., Vol. 6, pp.420–433, 1976.) Hummel and Zucker introduced a standard theory ofconsistency in labeling problems that established a connection withmathematical optimization techniques. (R. A. Hummel and S. W. Zucker,‘On the foundation of relaxation labeling processes’, IEEE Trans.Pattern Anal. Machine Intell., Vol. 5, pp 267–287, 1983.) Much of theliterature is concerned with application of standard techniques andrecently Chen and Luh (1995) and Pelillo (1997) have published papersthat cover extensions to the theoretical basis established by Hummel andZucker. (Qin Chen and J. Y. S. Luh, ‘Relaxation labeling algorithm forinformation integration and its convergence’, Pattern Recognition, Vol.28, No. 11, pp. 1705–1722, 1995; M. Pelillo, ‘On the dynamics ofnonlinear relaxation labeling processes’, Journal of MathematicalImaging and Vision 7, pp. 309–323, 1997, Kluver Academic Publishers.)Some of the early work is summarized by Kittler and Illingworth. [J.Kittler and J. Illingworth, ‘Relaxation labeling algorithms—a review’,Image Vision Comp., Vol. 3, pp. 206–215, 1985.]

The relaxation labeling process involves a number of entities. There isa set of N objects b_(i)εB to be labeled with a set λ={1,2,3, . . . ,M}of M labels. In the notation, Roman letters represent indices forobjects and Greek letters the indices for labels. Weighted labelingassignments P_(i)(λ) provide a certainty measure for assigning orassociating labels λ to object b_(i). Because P_(i)(λ) satisfies thefollowing probability axioms0≦P _(i)(λ)≦1.0,all i,λ,  Eqn. 2.1

$\begin{matrix}{{{\sum\limits_{\lambda = 1}^{M}{P_{i}(\lambda)}} = 1.0},\mspace{14mu}{i = 1},\ldots\mspace{14mu},N,} & {{Eqn}.\mspace{14mu} 2.2}\end{matrix}$

it may be interpreted as a probability. Different authors refer toP_(i)(λ) differently. For example, it is sometimes a certainty vector ora weighting vector. In this paper, P_(i)(λ) is called the assignmentvector. The concatenation of assignment vectors for all objects in B isdenoted by P≡{right arrow over (P)} (M×N elements).

The space of assignment vectors consistent with 2.1 and 2.2 is denotedby IK. Every vertex of IK represents an unambiguous labeling assignmentand this subspace is denoted by IK *. A labeling P, which is an interiorpoint of IK, is called strictly ambiguous. The iteration index K inP^(K) _(I)(λ) is generally suppressed in order to simplify notation. Theinitial distributions may be taken to be a uniform distribution havingmaximum entropy; i.e.,

$\begin{matrix}{{{P_{i}^{K = 0}(\lambda)} = \frac{1}{M}},} & {{Eqn}.\mspace{14mu} 2.3}\end{matrix}$

that is suitable if there is no specific information available formaking the initial assignment of labels.

A compatibility matrix R_(ij)(α,β) is an array of compatibilitycoefficients or consistency constraints, which are defined in terms ofsensor reports and other available intelligence. The compatibilitycoefficients define the relative influence of object j having label β onthe assignment of object i with label α. For example,R_(ij)(α,β)>R_(ij)(α,γ) implies that label β is more compatible thanlabel γ for object j whenever label α is assigned to object i. The goalof the labeling process is to assign labels to each object (i.e., solvefor the assignment vectors) in a manner that is consistent with respectto the compatibility constraints. A dependent matrix called the supportfunction is defined by

$\begin{matrix}{{S_{i}(\lambda)} = {\sum\limits_{j = 1}^{N}{\sum\limits_{\lambda^{\prime} = 1}^{M}{{R_{ij}\left( {\lambda,\lambda^{\prime}} \right)}{{P_{j}\left( \lambda^{\prime} \right)}.}}}}} & {{Eqn}.\mspace{14mu} 2.4}\end{matrix}$

Defining the compatibility matrix is a modeling task separate from thetask of finding a “consistent” solution for the assignment of labels.The RLA provides one well studied mathematical method for solving forthe assignment vector and is based on a non-linear update rule.Historically, we have two different versions of the update rule:

$\begin{matrix}{{{P^{K + 1}{i(\lambda)}} = \frac{{P_{i}^{K}(\lambda)} \cdot \left\{ {1 + {S_{i}^{K}(\lambda)}} \right\}}{1 + {{\hat{S}}^{K}i}}},\text{[Rosenfeld,~~Hummel~~and~~Zucker~~1971]}} & {{Eqn}.\mspace{14mu} 2.5} \\{{{{where}\mspace{14mu}{\hat{S}}^{K}i} = {\sum\limits_{\lambda = 1}^{M}{{P_{i}^{K}(\lambda)} \cdot {S_{i}(\lambda)}}}},\mspace{14mu}{and}} & {{Eqn}.\mspace{14mu} 2.6} \\{{P^{K + 1}{i(\lambda)}} = {\frac{{P_{i}^{K}(\lambda)} \cdot {S_{i}^{K}(\lambda)}}{{\hat{S}}_{i}^{K}(\lambda)}.}} & {{Eqn}.\mspace{14mu} 2.7}\end{matrix}$

[S. W. Zucker, Y. G. Leclerc, and J. L. Mohammed, ‘Continuous relaxationand local maxima selection: conditions for equivalence’, IEEE Trans.Pattern Anal. Mach. Intell., PAMI-3, pp. 117–137, March 1981.]

Note that the normalization insures that the probability axioms 2.1 and2.2 are satisfied.

The following is a brief review of material relating to theHummel-Zucker consistency theory that is covered more thoroughly inPelillo's paper. The relaxation algorithm may be viewed as a continuousmapping of the assignment space onto itself and uses the notationP ^(K+1)=

(P ^(K)),K≧0,PεIK,  Eqn. 2.8

that continues the process until it reaches a fixed or equilibriumpoint, whereP ^(K)=

(P ^(K)),for some K,  Eqn. 2.9

and S_(i)(λ)=c_(i) whenever P_(i)(λ)>0, i=1,2, . . . ,N,λεL for somenonnegative constants c_(i).

A weighted labeling assignment PεIK is said to be consistent if thefollowing N inequalities are satisfied,

$\begin{matrix}{{{\sum\limits_{\lambda = 1}^{M}{{P_{i}(\lambda)} \cdot {S_{i}\left( {\lambda;\overset{\rightarrow}{P}} \right)}}} \geq {\sum\limits_{\lambda = 1}^{M}{{V_{i}(\lambda)} \cdot {S_{i}\left( {\lambda;\overset{\rightarrow}{P}} \right)}}}}{i = 1},2,\ldots\mspace{14mu},\;{N\mspace{14mu}{for}\mspace{14mu}{all}\mspace{14mu} V\mspace{14mu}{in}\mspace{14mu}{{IK}.}}} & {{Eqn}.\mspace{14mu} 2.10}\end{matrix}$

If strict inequalities hold in 2.10, then P is said to be strictlyconsistent.

Theorem [3.1 Pelillo]: A labeling PεIK is consistent if and only if forall i=1,2, . . . ,N the following conditions hold:S _(i)(λ)=c _(i),whenever P _(i)(λ)>0  1)S _(i)(λ)≦c _(i),whenever P _(i)(λ)=0  2)

for some nonnegative constants c_(i).

Corollary [3.2 Pelillo]: Let PεIK be consistent. Then P is a fixed pointfor the nonlinear relaxation operator ℑ. Moreover, if P is strictlyambiguous the converse also holds.

Hummel and Zucker introduced the average local consistency, defined as

$\begin{matrix}{{A\left( \overset{\rightarrow}{P} \right)} = {\sum\limits_{i = 1}^{N}{\sum\limits_{\lambda = 1}^{M}{{P_{i}(\lambda)} \cdot {S_{i}(\lambda)}}}}} & {{Eqn}.\mspace{14mu} 2.11}\end{matrix}$

and proved the following theorem:

Theorem [Hummel-Zucker] If the compatibility matrix R is symmetric;i.e., R_(ij)(α,β)=R_(ij)(α,β), for all i,j,α,β, then any local maximumPεIK of A is consistent.

Note that the inverse of the theorem need not be true.

The requirement for symmetry allows one to relate the gradient of theaverage local consistency to the support function

$\begin{matrix}{\frac{\partial{A\left( \overset{\rightarrow}{P} \right)}}{\partial{p_{i}(\lambda)}} = {{2{\sum\limits_{j = 1}^{N}{\sum\limits_{\beta = 1}^{M}{{R_{ij}\left( {\lambda,\beta} \right)} \cdot {P_{j}(\beta)}}}}} = {2{{S_{i}(\lambda)}.}}}} & {{Eqn}.\mspace{14mu} 2.12}\end{matrix}$

The compatibility matrix need not be restricted to interactions betweenpairs of labeled objects, but could entail higher order relaxationschemes such as the third order compatibility relations that follow

$\begin{matrix}{{S_{i}(\lambda)} = {\sum\limits_{j = 1}^{N}{\sum\limits_{\beta = 1}^{M}{\sum\limits_{k = 1}^{N}{\sum\limits_{\gamma = 1}^{M}{{R_{ijk}\left( {\alpha,\beta,\gamma} \right)}{P_{j}(\beta)}{{P_{k}(\gamma)}.}}}}}}} & {{Eqn}.\mspace{14mu} 2.13}\end{matrix}$

If the compatibility matrix satisfies the third order symmetryrestrictionR _(ijk)(α,β,γ)+R _(kij)(γ,α,β)+R _(jki)(β,γ,α)=3R _(ijk)(α,β,γ),  Eqn.2.14

then the gradient satisfies

$\begin{matrix}{\frac{\partial{A\left( \overset{\rightarrow}{P} \right)}}{\partial{p_{i}(\lambda)}} = {3{{S_{i}(\lambda)}.}}} & {{Eqn}.\mspace{14mu} 2.15}\end{matrix}$

Although an asymmetric compatibility matrix provides no objectivefunction such as A({right arrow over (P)}) to be maximized, therelaxation labeling algorithm is still useable. The application ofrelaxation algorithms to arbitrary compatibility matrices is thoroughlydiscussed in Pelillo (1997).

Another influence on the RLA formalism used in A-Track comes from Chenand Luh (1995) which presents new formats for the compatibility matrixand a new updating algorithm to replace standard updating equations 2.5and 2.7. Their format involves the use of a compatibility block matrixof the form

$\begin{matrix}{C = \begin{bmatrix}{\underset{\_}{C}}_{11} & {\underset{\_}{C}}_{12} & \ldots & {\underset{\_}{C}}_{1N} \\{\underset{\_}{C}}_{21} & \ldots & \ldots & \ldots \\\ldots & \ldots & \ldots & \ldots \\{\underset{\_}{C}}_{N1} & {\underset{\_}{C}}_{N2} & \ldots & {\underset{\_}{C}}_{NN}\end{bmatrix}} & {{Eqn}.\mspace{14mu} 2.16}\end{matrix}$

where C _(ij)εIR^(M×M) and CεIR^(MN×MN). Matrices are denoted byunderlining or bold script. A neighborhood weight matrix is defined by

$\begin{matrix}{\underset{\_}{d} = \begin{bmatrix}d_{11} & d_{12} & \ldots & d_{1M} \\d_{21} & \ldots & \ldots & \ldots \\\ldots & \ldots & \ldots & \ldots \\d_{M1} & d_{M2} & \ldots & d_{MM}\end{bmatrix}} & {{Eqn}.\mspace{14mu} 2.17}\end{matrix}$

where

$\begin{matrix}{{{0 \leq d_{ij} \leq {1.0\mspace{14mu}{and}\mspace{14mu}{\sum\limits_{j = 1}^{M}d_{ij}}}} = {1.0\mspace{14mu}{for}\mspace{14mu} i}},{j = 1},2,\ldots\mspace{14mu},{M.}} & {{Eqn}.\mspace{14mu} 2.18}\end{matrix}$

The support matrix is given by

$\begin{matrix}{{\overset{\rightarrow}{S}}_{i} = {\sum\limits_{j = 1}^{M}{d_{ij} \cdot {{{\underset{\_}{C}}_{ij} \circ {\overset{\rightharpoonup}{P}}_{j}}.}}}} & {{Eqn}.\mspace{14mu} 2.19}\end{matrix}$

The new updating formula isP _(i) ^(K+1)(λ)=P _(i) ^(K)(λ){1+S _(i) ^(K)(λ)−Ŝ _(i) ^(K)}  Eqn. 2.20

that is based on the following two conditions on the compatibilitycoefficients:

$\begin{matrix}{{{\sum\limits_{\alpha = 1}^{M}{C_{ij}\left( {\alpha,\beta} \right)}} = {1.0\mspace{14mu}{where}}}{i,{j = 1},\ldots\mspace{14mu},{N\mspace{14mu}{and}\mspace{14mu}\alpha},{\beta = 1},\ldots\mspace{14mu},M,}} & {{{Eqn}.\mspace{14mu} 2.21}a}\end{matrix}$and 0≦C _(ij)(α,β)≦1.0.  Eqn. 2.21b

With conditions 2.21, one can verify that

$\begin{matrix}{{\sum\limits_{\alpha = 1}^{M}{S_{i}(\alpha)}} = {{1.0\mspace{14mu}{and}\mspace{14mu} 0} \leq {S_{i}(\alpha)} \leq 1.0}} & {{Eqn}.\mspace{14mu} 2.22}\end{matrix}$

and0≦Ŝ_(i)≦1.0.  Eqn. 2.23

Recall equation 2.6,

${\hat{S}}_{i}^{K} = {\sum\limits_{\lambda = 1}^{M}{{P_{i}^{K}(\lambda)} \cdot {{S_{i}^{K}(\lambda)}.}}}$It can be shown also that P^(K+1) satisfies the probability axioms,equations 2.1 and 2.2, if P^(K) does. The normalization employed in theother updating formulas (equations 2.5 and 2.7) to insure that P^(K+1)satisfies the probability axioms is no longer required provided theinitial vector P⁰ satisfies the axioms.

To summarize the updating equations with slightly different notation,letdP ^(K) _(i)(λ)=P ^(K+1) _(i)(λ)−P ^(K) _(i)(λ)  Eqn. 2.24

anddS ^(K) _(i)(λ)=S ^(K) _(i)(λ)−Ŝ ^(K) _(i).  Eqn. 2.25

Then we havedP ^(K) _(i)(λ)=P ^(K) _(i)(λ)·ds _(i) ^(K)(λ),[Chen and Luh 1995]  Eqn.2.26

$\begin{matrix}{{{{dP}_{i}^{K}(\lambda)} = \frac{{P_{i}^{K}(\lambda)} \cdot {{dS}_{i}^{K}(\lambda)}}{\left( {1 + {\hat{S}}_{i}^{K}} \right)}},\mspace{14mu}\text{[Rosenfeld,~~et~~al.~~1978]}} & {{Eqn}.\mspace{14mu} 2.27} \\{{{dP}_{i}^{K}(\lambda)} = {\frac{{P_{i}^{K}(\lambda)} \cdot {{dS}_{i}^{K}(\lambda)}}{{\hat{S}}_{i}^{K}}.\mspace{20mu}\text{[Zucker,~~et~~al.~~1981]}}} & {{Eqn}.\mspace{14mu} 2.28}\end{matrix}$

In the next section, we adopt a modified version of the Chen and Luhupdate algorithm that differs from the compatibility matrix constraintsexpressed in equations 2.21a and 2.21b.

3. Modified Relaxation Labeling Algorithm (RLA) for A-Track

Generally known standard relaxation labeling algorithms may beadvantageously modified before they can be applied to A-Track 10. InA-Track 10, the objects are tracks denoted by lower case Roman letters.A track conceptually links a group of reports associated with separateframes. Therefore, dual indices are used to specify a specific sensorreport: a frame index (capital Greek letter) and a report index (lowercase Greek letter) within a frame. The assignment vector P(i,Ω,α) forthe i^(th) track in the Ω-frame gives the relative assignment weightsfor the α-indexed reports in the Ω-frame.

The A-Track system 10 uses both second and third order compatibilitymatrices:(2^(nd) order)R(i,Ω,α|j,Γ,β)  Eqn. 3.1(3^(rd) order)R(i,Ω,α|j,Γ,β|,k,Λ,γ)  Eqn. 3.2

We modified the notation of Section 2 because of the additional frameindex needed for labeling; for example, P(i,Ω,α) is shorthand forP_(Ω,i)(Ω,α) and R(i,Ω,α|j,Γ,β) for R_(Ω,i;Γ,j)(Ω, α;γ,γ). With a littlepractice, the modified notation becomes easier to read and write. Thinkof the frame index as shared between the track and report identifiers.

Support function definitions require extra summations over the frameindex:

$\begin{matrix}{{S\left( {i,\Omega,\alpha} \right)} = {{\sum\limits_{j = 1}^{N}{\sum\limits_{\Gamma = 0}^{F}{\sum\limits_{\beta = 0}^{M_{\Gamma}}{{R\left( {i,\Omega,\left. \alpha \middle| j \right.,\Gamma,\beta} \right)}{P\left( {j,\Gamma,\beta} \right)}}}}} \equiv {\sum\limits_{j,\Gamma,\beta}{{R\left( {i,\Omega,\left. \alpha \middle| j \right.,\Gamma,\beta} \right)}{P\left( {j,\Gamma,\beta} \right)}\mspace{14mu}{or}}}}} & {{Eqn}.\mspace{14mu} 3.3} \\{\mspace{79mu}{{S\left( {i,\Omega,\alpha} \right)} = {\sum\limits_{j,\Gamma,\beta}{d_{ij}D_{\Omega\;\Gamma}{R\left( {i,\Omega,\left. \alpha \middle| j \right.,\Gamma,\beta} \right)}{P\left( {j,\Gamma,\beta} \right)}}}}} & {{Eqn}.\mspace{14mu} 3.4}\end{matrix}$

where the latter version uses dual neighborhood weight matrices whichsatisfy

$\begin{matrix}{{{{0 \leq d_{ij} \leq {1.0\mspace{14mu}{and}\mspace{14mu}{\sum\limits_{j}^{N}d_{ij}}}} = {1.0\mspace{14mu}{for}\mspace{14mu} i}},{j = 1},2,\ldots\mspace{14mu},N}\;{and}} & {{Eqn}.\mspace{14mu} 3.5} \\{{{0 \leq D_{\Omega\;\Gamma} \leq {1.0\mspace{14mu}{and}\mspace{14mu}{\sum\limits_{\Gamma = 0}^{F}D_{\Omega\;\Gamma}}}} = {1.0\mspace{14mu}{for}}}\text{}{\Omega,{\Gamma = 0},1,2,\ldots\mspace{14mu},{F\;.}}} & {{Eqn}.\mspace{14mu} 3.6}\end{matrix}$

Each frame has M_(Γ) sensor reports denoted by the indices {1,2, . . .,MΓ} and a zero label that denotes a missing sensor report in a givenframe. A-Track employs two sliding time windows, each of which containsa fixed number of frames. One is an active window in which theassignment vectors are unknown variables. Note that in Equations 3.3 and3.4, the support function is defined in terms of the F+1 frames in theactive window. The second sliding window is associated with historicframes that have most recently passed out of the active window.Assignment vectors in the historic frames are frozen and no longersubject to variation. The incorporation of historic frames in theA-Track model is addressed in Section 4.

In A-Track, the number of objects (tracks) is an unknown and varies overtime as old tracks disappear and new tracks are created. In the standardRLA applications, the object and label sets are generally fixed sizes.Here the number of sensor reports associated with each frame is fixedbut varies for different frames. The creation and initialization of newtracks are discussed in Section 7. It can be noted here that thecreation of a new track does not alter existing coefficients of thecompatibility matrices; however, it does require the computation of newcoefficients for the new track. In order to simplify this presentationit is assumed that there is at most one report per track in each frameand defer further discussion of the multiple report case to Section 4.

The remainder of this section is devoted to developing the formalism fora new RLA that utilizes Chen and Luh's updating formula in equation 2.20or 2.26. In the Chen-Luh algorithm, the compatibility matrix mustsatisfy two constraints. Equation 2.21b requires the compatibilitymatrix to have a positive dynamic range in the interval [0,1] thatinsures the invariance of the probability axioms for P(i,Ω,α).

Constraint 2.21a, which normalizes the compatibility matrix, isnecessary to insure that the support function also satisfies theprobability axioms. The disadvantage of this normalization is that itremoves any symmetry in the compatibility matrix. The new RLA formalismdeveloped for A-Track preserves the simplified update formula withoutsacrificing symmetry in the compatibility matrix and also facilitatesthe use of negative coefficients and support functions that depend on amultiplicity of compatibility matrices with mixed orders.

To combine a set of compatibility matrices we take a weighted average ofthe support functions defined with respect to each compatibility matrix;i.e.,

$\begin{matrix}{\mspace{79mu}{{S\left( {i,\Omega,\alpha} \right)} = {\sum\limits_{Z}{{W_{Z} \cdot {S_{Z}\left( {i,\Omega,\alpha} \right)}}\mspace{14mu}{where}}}}} & {{Eqn}.\mspace{14mu} 3.7} \\{{S_{Z}\left( {i,\Omega,\alpha} \right)} = {\sum\limits_{j,\Gamma,{\beta\ldots}}{{{R_{Z}\left( {i,\Omega,\left. \alpha \middle| i \right.,\Gamma,{\beta\mspace{14mu}\ldots}} \right)} \cdot {P\left( {j,\Gamma,\beta} \right)}}\mspace{14mu}\ldots\mspace{14mu}{and}}}} & {{Eqn}.\mspace{14mu} 3.8} \\{\mspace{79mu}{{\sum\limits_{Z}W_{Z}} = {1.0.}}} & {{Eqn}.\mspace{14mu} 3.9}\end{matrix}$

The Z index denotes elements in the sets of support functions andcompatibility matrices. The average local consistency defined inequation 2.22 is re-defined as

$\begin{matrix}{{{{A\left( \overset{\rightarrow}{P} \right)} =}\quad}{\quad{\sum\limits_{Z}{\sum{\sum{W_{Z} \cdot {R_{Z}\left( {i,\Omega,\left. \alpha \middle| j \right.,\Gamma,{\beta\mspace{14mu}\ldots}} \right)} \cdot {P\left( {i,\Omega,\alpha} \right)} \cdot}}}}\quad}{P\left( {j,\Gamma,\beta} \right)}\mspace{14mu}\ldots} & {{Eqn}.\mspace{14mu} 3.10}\end{matrix}$

and is the A-Track objective function to be maximized. If all thecompatibility matrices are symmetric, then the gradient has the form

$\begin{matrix}{{\frac{\partial{A\left( \overset{\rightarrow}{P} \right)}}{\partial{P\left( {i,\Omega,\alpha} \right)}} = {\sum\limits_{Z}{W_{Z} \cdot O_{Z} \cdot {S_{Z}\left( {i,\Omega,\alpha} \right)}}}},} & {{Eqn}.\mspace{14mu} 3.11}\end{matrix}$

where O_(Z) equals the ‘effective’ order of the compatibility matricesR_(Z). The effective order is defined as the number of variable Ps inthe objective function (Equation 3.10) which are in the active slidingwindow. In this paper, O_(Z) is a positive integer between 1 and 3.Recall that the Ps associated with historic frames are considered fixedparameters or inactive assignment vectors.

The new updating formalism is contained in the following three theorems.The conditions in the first theorem will be assumed in the enhanced RLAused in A-Track.

Theorem 3.1 If given the following three conditions−1.0≦L ^(Z) _(MIN) ≦R _(Z)(i,Ω,α|j,Γ,β . . . )≦L ^(Z) _(MAX)≦1.0 for allZ,i,Ω,α,j,Γ,β, . . . ,  Eqn. 3.12L ^(Z) _(MAX) −L _(MIN)≦1.0,  Eqn. 3.13

$\begin{matrix}{{0 \leq {P\left( {i,\Omega,\alpha} \right)} \leq {1.0\mspace{14mu}{and}}}\;{{{\sum\limits_{\alpha}{P\left( {i,\Omega,\alpha} \right)}} = {1.0\mspace{14mu}\left( {{probability}\mspace{14mu}{axioms}} \right)}},}} & {{Eqn}.\mspace{14mu} 3.14}\end{matrix}$

then the following results are self evident:L ^(Z) _(MIN) ≦S _(Z)(i,Ω,α)≦L ^(Z) _(MAX)  Eqn. 3.15

$\begin{matrix}{L_{MIN}^{Z} \leq {{\hat{S}}_{Z}\left( {i,\Omega,\alpha} \right)} \equiv {\sum\limits_{\alpha}{{S_{Z}\left( {i,\Omega,\alpha} \right)} \cdot {P\left( {i,\Omega,\alpha} \right)}}} \leq L_{MAX}^{Z}} & {{Eqn}.\mspace{14mu} 3.16}\end{matrix}$−1.0≦dS _(Z)(i,106 ,α)≡S _(Z)(i,Ω,α)−Ŝ _(Z)(i,Ω,α)≦1.0  Eqn. 3.17

$\begin{matrix}{{- 1.0} \leq {\mathbb{d}{S\left( {i,\Omega,\alpha} \right)}} \equiv {\sum\limits_{Z}{W_{Z} \cdot {\mathbb{d}{S_{Z}\left( {i,\Omega,\alpha} \right)}}}} \leq 1.0} & {{Eqn}.\mspace{14mu} 3.18} \\{where} & \; \\{{\sum\limits_{Z}W_{Z}} = {{1.0\mspace{14mu}{and}\mspace{14mu} 0} \leq W_{Z} \leq {1.0.}}} & {{Eqn}.\mspace{14mu} 3.19}\end{matrix}$

Proof of these assertions only requires substitution of the limits intothe basic definitions. The shorthand notation of this section can beeasily converted into the standard notation describe in Section 2 foruse in other problems not requiring multiple labeling indices. Similarremarks hold for the remaining theorems in this section.

The second theorem modifies the Chen-Luh RLA by replacing thenormalization conditions on the compatibility matrix given by equation2.21b with conditions 3.12 and 3.13 in Theorem 1 and completelyeliminating their normalization constraints on the compatibilitycoefficients expressed in equation 2.21a.

Theorem 3.2: If the three conditions in Theorem 3.1 are valid anddP ^(K)(i,Ω,α)=P ^(K)(i,Ω,α)·dS ^(K)(i,Ω,α)for all i,Ω,α  Eqn. 3.20

where

K=iteration index,

$\begin{matrix}{{{S^{K}\left( {i,\Omega,\alpha} \right)} \equiv {\sum\limits_{Z}{W_{Z} \cdot {S_{Z}^{K}\left( {i,\Omega,\alpha} \right)}}}},} & {{Eqn}.\mspace{14mu} 3.21} \\{{0 \leq W_{Z} \leq 1},\mspace{14mu}{{\sum\limits_{Z}W_{Z}} = 1.0},} & {{Eqn}.\mspace{14mu} 3.22}\end{matrix}$dS ^(K)(i,Ω,α)≡S ^(K)(i,Ω,α)−Ŝ ^(K)(i,Ω)for all i,Ω,α,  Eqn. 3.23

$\begin{matrix}{{{{\hat{S}}^{K}\left( {i,\Omega} \right)} \equiv {\sum\limits_{\alpha}{\sum\limits_{Z}{W_{Z} \cdot {S_{Z}^{K}\left( {i,\Omega,\alpha} \right)} \cdot {P^{K}\left( {i,\Omega,\alpha} \right)}}}} \equiv {\sum\limits_{\alpha}{{S^{K}\left( {i,\Omega,\alpha} \right)} \cdot {P^{K}\left( {i,\Omega,\alpha} \right)}}}},} & {{Eqn}.\mspace{14mu} 3.24}\end{matrix}$

then the simplified update rule 3.20 preserves the probability axioms;i.e.,0≦P ^(K+1)(i,Ω,α)≡P ^(K) +dP ^(K)(i,Ω,α)≦1.0  Eqn. 3.25

and

$\begin{matrix}{{\sum\limits_{\alpha}{P^{K + 1}\left( {i,\Omega,\alpha} \right)}} = {1.0.}} & {{Eqn}.\mspace{14mu} 3.26}\end{matrix}$

Proof

By definition of Ŝ(i,Ω), there is the following identity,

$\begin{matrix}{{\sum\limits_{\alpha}{\mathbb{d}{P^{k}\left( {i,\Omega,\alpha} \right)}}} \equiv {\sum\limits_{\alpha}{{P^{K}\left( {i,\Omega,\alpha} \right)} \cdot {\mathbb{d}{S^{K}\left( {i,\Omega,\alpha} \right)}}}} \equiv 0} & {{Eqn}.\mspace{14mu} 3.27}\end{matrix}$

where

${{\sum\limits_{\alpha}{\mathbb{d}{P^{K}\left( {i,\Omega,\alpha} \right)}}} \equiv {{\sum\limits_{\alpha}{{P^{K}\left( {i,\Omega,\alpha} \right)} \cdot {S^{K}\left( {i,\Omega,\alpha} \right)}}} - {{{\hat{S}}^{K}\left( {i,\Omega} \right)}{\sum\limits_{\alpha}{P^{K}\left( {i,\Omega,\alpha} \right)}}}}} = {{{{\hat{S}}^{K}\left( {i,\Omega} \right)} - {{\hat{S}}^{K}\left( {i,\Omega} \right)}} \equiv 0.}$

Equations 3.27 and 3.14 establish equation 3.26. Substitution ofequation 3.18 (Theorem 1),−1.0≦dS ^(K)(i,Ω,α)≡S ^(K)(i,Ω,α)−Ŝ ^(K)(i,Ω)≦1.0 for all i,Ω,α,  Eqn.3.18

into the update rule 3.20 yields−P^(K)(i,Ω,α)≦dP ^(K)(i,Ω,α)≦P ^(K)(i,Ω,α).  Eqn. 3.28

Hence,P ^(K+1)(i,Ω,α)=P ^(K)(i,Ω,α)+dP ^(K)(i,Ω,α)≧0.  Eqn. 3.29

The upper limit of equation 3.25 follows immediately from equations 3.26and 3.29. QED

The third theorem gives a specific form for the simplified updatealgorithm that depends on the gradient of the average local consistency.

Theorem 3.3 If all compatibility matrices are symmetric and satisfy theconditions of Theorem 1 and

$\begin{matrix}{{S\left( {i,\Omega,\alpha} \right)} = {\sum\limits_{Z}{W_{Z} \cdot {S_{Z}\left( {i,\Omega,\alpha} \right)}}}} & {{Eqn}.\mspace{14mu} 3.30}\end{matrix}$

and the update algorithmdP(i,Ω,α)=P(i,Ω,α)·dG(i,Ω,α)  Eqn. 3.31

is dependent on G(i,Ω,α) which is proportional to the normalizedgradient of the average local consistency

$\begin{matrix}\left. {{A\left( \overset{\rightarrow}{P} \right)} = {{\sum\limits_{\alpha}{{S\left( {i,\Omega,\alpha} \right)} \cdot {P\left( {i,\Omega,\alpha} \right)}}} = {\sum\limits_{Z}{\sum\limits_{\alpha}{W_{Z} \cdot {S_{Z}\left( {i,\Omega,\alpha} \right)} \cdot {P\left( {i,\Omega,\alpha} \right)}}}}}} \right) & {{Eqn}.\mspace{14mu} 3.32} \\{where} & \; \\{{\frac{\partial{A\left( \overset{\rightarrow}{P} \right)}}{\partial{P\left( {i,\Omega,\alpha} \right)}} = {{\sum\limits_{Z}{W_{Z} \cdot O_{Z} \cdot {S_{Z}\left( {i,\Omega,\alpha} \right)}}} = {{G\left( {i,\Omega,\alpha} \right)} \cdot {\sum\limits_{Z}{W_{Z} \cdot O_{Z}}}}}},} & {{Eqn}.\mspace{14mu} 3.33} \\{{{\mathbb{d}{G\left( {i,\Omega,\alpha} \right)}} = {{{\sum\limits_{Z}{V_{Z} \cdot {S_{Z}\left( {i,\Omega,\alpha} \right)}}} - {\sum\limits_{\alpha}{{P\left( {i,\Omega,\alpha} \right)} \cdot {G\left( {i,\Omega,\alpha} \right)}}}} = {{G\left( {i,\Omega,\alpha} \right)} - {\hat{G}\left( {i,\Omega} \right)}}}},} & {{Eqn}.\mspace{14mu} 3.34} \\{and} & \; \\{{V_{Z} = \frac{W_{Z} \cdot O_{Z}}{\sum\limits_{Z}{W_{Z} \cdot O_{Z}}}},} & {{Eqn}.\mspace{14mu} 3.35}\end{matrix}$

then the derivative of A({right arrow over (P)}) is positive or zero:

$\begin{matrix}{{\mathbb{d}{A\left( \overset{\rightharpoonup}{P} \right)}} = {{\sum\limits_{i,\Omega,\alpha}{\frac{\partial{A\left( \overset{\rightharpoonup}{P} \right)}}{\partial{P\left( {i,\Omega,\alpha} \right)}} \cdot {\mathbb{d}{P\left( {i,\Omega,\alpha} \right)}}}} \geq {0.0.}}} & {{Eqn}.\mspace{14mu} 3.36}\end{matrix}$

Proof

Substitution of equation 3.33 into equation 3.36 yields

$\begin{matrix}{{{\mathbb{d}A} \propto {\sum\limits_{i,\Omega,\alpha}{{G\left( {i,\Omega,\alpha} \right)} \cdot {\mathbb{d}{P\left( {i,\Omega,\alpha} \right)}}}}} = {\sum\limits_{i,{\Omega\alpha}}{{G\left( {i,\Omega,\alpha} \right)} \cdot {P\left( {i,\Omega,\alpha} \right)} \cdot {{\mathbb{d}{G\left( {i,\Omega,\alpha} \right)}}.}}}} & {{Eqn}.\mspace{14mu} 3.37}\end{matrix}$

Substitution of G(i,Ω,α)≡Ĝ(i,Ω)+dG(i,Ω,α) into equation 3.37 yields

$\begin{matrix}{{\mathbb{d}A} \propto {{{\hat{G}\left( {i,\Omega} \right)} \cdot {\sum\limits_{i,{\Omega\alpha}}{{P\left( {i,\Omega,\alpha} \right)} \cdot {\mathbb{d}{G\left( {i,\Omega,\alpha} \right)}}}}} + {\sum\limits_{i,\Omega,\alpha}{{\mathbb{d}{G\left( {i,\Omega,\alpha} \right)}} \cdot {P\left( {i,\Omega,\alpha} \right)} \cdot {\mathbb{d}{G\left( {i,\Omega,\alpha} \right)}}}}}} & {{Eqn}.\mspace{14mu} 3.38}\end{matrix}$

Use of the identity,

$\begin{matrix}{{\sum\limits_{\alpha}{{\mathbb{d}{G\left( {i,\Omega,\alpha} \right)}} \cdot {P\left( {i,\Omega,\alpha} \right)}}} \equiv 0.0} & \;\end{matrix}$(see equation 3.27) reduces equation 3.38 to

$\begin{matrix}{{{\mathbb{d}A} \propto {\sum\limits_{i,{\Omega\alpha}}{{P\left( {i,\Omega,\alpha} \right)} \cdot \left( {\mathbb{d}{G\left( {i,\Omega,\alpha} \right)}} \right)^{2}}} \geq 0.0},} & {{Eqn}.\mspace{14mu} 3.39}\end{matrix}$

which is obviously non-negative. QED

4. Basic Compatibility Matrix

4.1 Introduction

Having described the relaxation labeling process, the compatibilitymatrices are now defined. The objective here is to exploit all physicalmeasurements and other available and pertinent information such asdigital map data. The RLA formalism allows one to weight various typesof data and to form relative likelihoods or constraints expressed in thecompatibility matrix. Often in image processing applications, onecompatibility matrix is sufficient to capture the fundamentalconstraints used in the RLA. However, the HRR A-Track implementationdescribed in the following sections requires more than one compatibilitymatrix to accomplish this task.

This section defines the basic second order compatibility matrix,R(i,.Ω.α|j,Γ,β), which is the most complicated of all the compatibilitymatrices used in the A-Track HRR model. It is also used in thedefinition of another third order compatibility matrices in Section 8.This description, while elaborate in scope, is not meant to be theexclusive definition for the basic compatibility matrix, but to serve asan illustrative approach for the A-Track model.

The basic compatibility matrix must satisfy the following constraints:R(i,.Ω.α|j,Γ,β)=0 if i≠j, Same-track  Eqn. 4.10≦R(j,.Ω.α⊕j,Γ,β)≦1.0, Normalized  Eqn. 4.2R(i,.Ω.α|i,Γ,β)=R(i,Γ,β|i,.Ω.α), Symmetric  Eqn. 4.3

whereΩ ₀≦Ω,Γ≦Ω₀ +F,  Eqn. 4.4

F+1=number of frames in the active window,

Ω₀=index of oldest frame in the active window.

To reiterate, the basic matrix R(j,.Ω.α|j,Γ,β) represents the relativecompatibility or influence of association label β in frame Γ with trackj when label α in frame Ω is associated with the same track index j.Also, as discussed in the previous section, the symmetry constraintinsures that the gradient of average local consistency function ofR(i,.Ω.α|j,Γ,β) will be proportional to the active support function

$\begin{matrix}{{{V\left( {i,\Omega,\alpha} \right)} = {\sum\limits_{j}{\delta_{ij}{\sum\limits_{\Gamma = \Omega_{0}}^{\Omega_{0} + F}{D_{\Omega\Gamma}{\sum\limits_{\beta = 0}^{M_{\Gamma}}{{R\left( {i,\Omega,\left. \alpha \middle| j \right.,\Gamma,\beta} \right)} \cdot {P\left( {j,\Gamma,\beta} \right)}}}}}}}},} & {{Eqn}.\mspace{14mu} 4.5}\end{matrix}$

where

${\sum\limits_{\Gamma = \Omega_{0}}^{\Omega_{0} + F}D_{\Omega\Gamma}} = 1.0$and δ_(ij) is the Dirac delta function.

The average local consistency is given by

$\begin{matrix}\begin{matrix}{J_{V} = {\sum\limits_{i,\Omega,\alpha}{{V\left( {i,\Omega,\alpha} \right)} \cdot \left( {P\left( {i,\Omega,\alpha} \right)} \right.}}} \\{= {\sum\limits_{i,\Omega,\alpha}{\sum\limits_{\Gamma = \Omega_{0}}^{\Omega_{0} + F}{D_{\Omega\Gamma}{\sum\limits_{\beta}{{R\left( {i,\Omega,{\alpha ❘i},\Gamma,\beta} \right)} \cdot {P\left( {i,\Gamma,\beta} \right)} \cdot {P\left( {i,\Omega,\alpha} \right)}}}}}}}\end{matrix} & {{Eqn}.\mspace{14mu} 4.6}\end{matrix}$

with the gradient

$\begin{matrix}{\frac{\partial J_{v}}{\partial{P\left( {i,\Omega,\alpha} \right)}} = {2 \cdot {{V\left( {i,\Omega,\alpha} \right)}.}}} & {{Eqn}.\mspace{14mu} 4.7}\end{matrix}$

In contrast to the active support function, a historic version of thesupport function, H(i,Ω,α) is introduced in which the basiccompatibility matrix links historic frames with an active frame Ω in theactive window:

$\begin{matrix}{{H\left( {i,\Omega,\alpha} \right)} = {\sum\limits_{j}{\delta_{ij}{\sum\limits_{\Gamma = {\Omega_{0} - N_{H}}}^{\Omega_{0} - 1}{D_{\Omega\Gamma}{\sum\limits_{\beta = 0}^{M_{L}}{{R\left( {i,\Omega,\left. \alpha \middle| i \right.,\Gamma,\beta} \right)} \cdot {P\left( {i,\Gamma,\beta} \right)}}}}}}}} & {{Eqn}.\mspace{14mu} 4.8}\end{matrix}$

where Ω₀≦Ω≦Ω₀F, Γ<Ω₀, and N_(H) is the number of inactive frames Becausethe P(i,Γ,β) terms in equation 4.8 are associated with historic frames,they are not active variables and are not subject to change as isP(i,Ω,α) which is in the active window. Therefore, the gradientcorresponding to equation 4.7 is

$\begin{matrix}{{\frac{\partial J_{H}}{\partial{P\left( {i,\Lambda,\lambda} \right)}} = {\frac{\partial{\sum\limits_{i,\Omega,\alpha}{{H\left( {i,\Omega,\alpha} \right)} \cdot {P\left( {i,\Omega,\alpha} \right)}}}}{\partial{P\left( {i,\Lambda,\lambda} \right)}} = {H\left( {i,\Lambda,\lambda} \right)}}},} & {{Eqn}.\mspace{14mu} 4.9}\end{matrix}$

which differs by a factor of 2 from the form of equation 4.7.

The A-Track system 10 depicted in FIG. 3 helps to illustrate how theframe inputs are combined with knowledge bases and converted into trackdescriptions. The derivation of the basic compatibility matrix isaccomplished in two steps; compute a linking matrix and then compute thecompatibility matrix. The method for computing the linking matrixdepends on the frame rate or time interval between frames. Section 4.2takes a look at generating the linking matrix under harsh conditionscharacterized by low frame rates, absence of ATR, and partial knowledgeof target state vectors. Section 4.3 introduces the ATR contributionsand discusses additional possibilities for constructing components ofthe linking matrix under more favorable conditions. A co-moving targetmodel is presented in Section 4.4. The conversion of the linking matrixinto the compatibility matrix is covered in Section 4.5 and severaladvanced implementation issues are discussed in Section 4.6.

4.2 Mobility Matrix as a Component of the Linking Matrix

When frames are separated by large time intervals there are many morepotential associations between targets than when the frames are closetogether. Furthermore, it more difficult or unfeasible to use theprincipal of continuity to estimate target state vectors for largeseparations. To generate the relative weights between targets indifferent frames it becomes necessary to exploit techniques thatintegrate basic ATR constraints derived from the HRR sensor signatures,physical constraints derived from analysis of digital maps, andbehavioral heuristics applied to moving vehicles both singly and ingroups.

Initially, we concentrate only on those components that deal with tworeports in different frames and defer the consideration of componentsinvolving missing reports to a later section. Reports are ‘observables’which provide physical bases for constructing constraints and relativelikelihood measures reflected in matrix form. For the time being, assumeR(i,.Ω.,α|i,Ω,β)=0 for all α and β,  Eqn. 4.10

which rules out multiple frame reports for the same track.

Satisfaction of this constraint could result from the sensor hardware orsoftware. Below, the special case where there are multiple targetreports in the same frame will be discussed separately as well the casewhere the report labels α and β are zero indicating missing reports in aframe.

Each report (Ω,α) contains the following minimal sensor data:

Target position r(Ω,α) or r=ρ·ē_(φ) or (ρ,φ) relative to the sensor

Radial velocity {dot over ( r(Ω,α)={dot over (ρ)}(Ω,α)·ē_(φ)(Ω,α)

High Range Resolution (HRR) radar target signature Ψ(Ω,α).

In addition to sensor and ATR data, assume that digital maps of theroadways and other environmental features are available. Obtaining thisinformation is part of the preparation for a more extensivecompatibility matrix necessary to improve the tracking accuracy. It isvery important to know if the targets are moving along roads or offroads and if off road are they moving along pathways that serve asuncharted roads. These uncharted pathways can be learned over time byobservation of traffic flow patterns. In general, a target achieves amuch greater distance if it moves along a road or a used pathway than itwould across open country with rough terrain. A full presentation ofthis mode of analysis is beyond the scope of this paper. In thisdiscussion, we will refer to the exploitation of this kind ofinformation as mobility analysis.

The first practical step in mobility analysis is to rough out apossibility matrix based on an ultimate velocity limit between any tworeports which identifies the subset of possible components to thecompatibility matrix. Let E(Ω,α|Λ,β)ε{0,1}, where E=1 if a trajectorybetween two reports is physically possible and E=0 otherwise. Next,consider the four possible cases with regards to given pair reports forwhich E=1:

F1=1.0 if on-road to on-road

F1=0.8 if off-road to off-road

F1=0.6 if off-road to on-road or if on-road to off-road

The value of factor F1 reflects the relative likelihoods one couldassign to each of the given cases. Using the principal of continuity, itwould generally be agreed that vehicles on roads would tend to stay onthe road, F1=1. Similarly, for the continuity off-road we could assign avalue of 0.8 where the value is smaller than 1.0 because being on theroad may be more desirable if they are available. The switch between onand off road may be given the smallest value because it is a break incontinuity. Note the symmetry between on and off road switching isnecessary to maintain the symmetry of the compatibility matrix.

Further use of the continuity principle would lead one to expect avehicle to travel in the same general direction rather than reversedirection. Both options may be possible, but if one were to assignrelative weights it would look something like:

F2=1.0 Target maintains general direction,

F2=0.7 Target reverses direction.

The exact factor assignments are not crucial in defining a compatibilitymatrix but the correct biases or constraints should be reflected in thevalues chosen. Those skilled in mobility analysis of a particularenvironment would assign more realistic values than those selected forillustration.

Finally, the target class would be a major influence in assigning valuesto the mobility factors under different conditions. For example, themobility of tanks off-road may far exceed cars whereas on a major roadthe car may have the greater mobility. The factors F1 and F2 may becombined differently for different target classes. In practice a moreextensive set of heuristics and analytic algorithms would be used toassemble a mobility matrix M (Ω,α|Λ,β|c), which links two separatereports and is customized to different target classes indexed by c=1,2,. . . ,Nc. If M (Ω,α|Λ,β|c)=0 for a given class c, then it is impossibleto link these two reports for that class. Therefore, after suchanalysis, if M (Ω,α|Λ,β|c)=0 for all values of c, then the correspondingcomponents of E(Ω,α|Λ,β) must be updated.

Computing of the mobility matrix is one of the most computationallyintense tasks in A-Track. This is because all possible trajectoriesbetween two given reports must be considered with respect to terrainfeatures and roads over possibly large surveillance areas. Moreelaborate analysis would require specialized computers to reduce thecomputation time to match the frame intervals. The military exploitationof map information, which has been underway since the early nineteenseventies, will provide the initial bases for the analytic softwareneeded to compute the mobility matrix.

Perhaps in the future, the use of mathematical morphology on specializedcomputers could be employed. Here the idea is to use generalizeddilation operators that have non-uniform expansion rates in differentdirections and which adapt to localized features in the digital terrainand road maps. By conditionally initiating a sequence of dilations fromthe location of a single report, this report can be linked bothtemporally and spatially to all other reports by the spreading wavefront of the dilation process.

4.3 ATR Inputs and the Linking Matrix

The Identification (ID) vector P(c|Ω,α) for report (Ω,α) is derived fromthe target signature by using an ATR algorithm and specifies therelative probability the report is classified by the class index, c=(0),1, 2, . . . , Nc. The identification vectors come in two flavorsdepending on whether or not it uses the zero class index P(0|Ω,α) tospecify the relative likelihood of an unknown class. Let the ID vectorsatisfy the probability axioms,

$\begin{matrix}{{\sum\limits_{c}^{N_{c}}{P\left( {\left. \underset{.}{c} \middle| \Omega \right.,\alpha} \right)}} = {{1.0\mspace{14mu}{and}\mspace{14mu} 0} \leq {P\left( {\left. c \middle| \Omega \right.,\alpha} \right)} \leq {1.0\mspace{20mu}{for}\mspace{14mu}{all}\mspace{14mu}{c'}{s.}}}} & \;\end{matrix}$

The shorthand notation P(Ω,α) for the ID vector is also used.

If more sensor data were available to augment the minimal data set, onecould enhance the derived ID vector. In the limit with enough additionalsensor data and efficient ATR algorithms, it is possible to not onlyspecify a target class but to label individual targets within the sameclass. Obviously this kind of identification would greatly simplify thetracking problem. In this discussion, it is assumed that target ID's areimperfect at best, including the case where a target may be in anunknown class. It should be noted that it is not necessary to havecomplete certainty about the target ID vector for it to be a usefulconstraint in the matching and linking of different reports in the RLA.

The application of ATR algorithms to the HRR signature will now befactored into the development of the basic linking matrix. Later, aderived ID vector P(c|Ω,α) will be used to convert the linking matrixinto a compatibility matrix. HRR signatures are very sensitive to thepose of the target relative to the sensor, so in general, it is notpossible to use the following correlation measure directly between thesignatures in different frames separated by a large time interval;

$\begin{matrix}{0 \leq {C\left( {\Omega,\left. \alpha \middle| \Gamma \right.,\beta} \right)} \equiv \frac{{\overset{\_}{\Psi}\left( {\Omega,\alpha} \right)} \cdot {\overset{\_}{\Psi}\left( {\Gamma,\beta} \right)}}{{{\overset{\_}{\Psi}\left( {\Omega,\alpha} \right)}} \cdot {{\overset{\_}{\Psi}\left( {\Gamma,\beta} \right)}}} \leq {1.0.}} & {{Eqn}.\mspace{14mu} 4.11}\end{matrix}$

ATR template matching algorithms utilizes a similar correlation measurebetween the measured signature and reference signatures for sets ofclasses and poses (Θ):

$\begin{matrix}{0 \leq {C\left( {\Omega,{\alpha;c},\Theta} \right)} \equiv \frac{{\Psi_{Meas}\left( {\Omega,\alpha} \right)} \cdot {\Psi_{Ref}\left( {\Omega,{\alpha;c},\Theta} \right)}}{{{\Psi_{Meas}\left( {\Omega,\alpha} \right)}} \cdot {{\Psi_{Ref}\left( {\Omega,{\alpha;c},\Theta} \right)}}} \leq 1.0} & {{Eqn}.\mspace{14mu} 4.12}\end{matrix}$

In the HRR template matching algorithms, the maximum correlation valueover all orientationsC(Ω,α;c,Θ _(MAX))=MAX_(Θ) [C(Ω,α;c,Θ)],  Eqn. 4.13

generates the ID vector and also estimates the target pose Λ_(MAX). Itis even possible to generate an ID vector with a zero (‘other’)component; see Mitchell and Westerkamp (1999) for one approach appliedto HRR signatures. [Mitchell, Richard A. and John J. Westerkamp, ‘RobustStatistical Feature Based Aircraft Identification’, IEEE Transactions onAerospace and Electronic Systems, Vol. 35, No. 3, 1999.] Knowing thetarget pose can be utilized to generate dynamic constraints and targetvelocity vectors. Recall that target orientation together with themeasured radial velocity component provides an estimate of theinstantaneous velocity vector. Additional pose information can be usedto help confirm or improve the target ID and to improve the definitionsof co-moving target model describe below.

The main purpose of this section is to demonstrate how ATR informationis incorporated into A-Track. Subsequently, several ATR factors will bedefined that could be use with the mobility matrix to formulate alinking matrix. The first factor is a dot product of two ID vectors,which is a very general construct. It has the advantage of beinginvariant over pose variation. One difficulty with matching twodifferent ID vectors is how to handle the unknown class component (c=0)since this component can represent anything different from one of thegiven classes (c>0). In this discussion, we assign only partial creditto the matching of the unknown ‘zero’ component. If one thinks of the IDvectors as a probability distribution, then probability associated withmatching two ID vectors is the customized dot product (•) between them:

$\begin{matrix}{{{{\overset{\rightharpoonup}{P}\left( {\Omega,\alpha} \right)}\bullet{\overset{\rightharpoonup}{P}\left( {\Lambda,\beta} \right)}} = {\sum\limits_{c = 0}^{N_{c}}{\left\{ {1 - {v \cdot {\delta(c)}}} \right\} \cdot {P\left( {\left. c \middle| \Omega \right.,\alpha} \right)} \cdot {P\left( {\left. c \middle| \Lambda \right.,\beta} \right)}}}},} & {{Eqn}.\mspace{14mu} 4.14}\end{matrix}$Prob{ID(Ω,α)=c and ID(Λ,β)=c}≡P(c|Ω,α)·P(cβΛ,β) for c>0,  Eqn. 4.15Prob{ID(Ω,α)=0 and ID(Λ,β)=0}≡{1−v}·P(0|Ω,α)·P(0|Λ,β) for c=0,  Eqn.4.16where δ(c) is the Dirac delta function (δ(c)=1 if c=0 and =0 otherwise)and 0<v<1. Combining the matching coefficients and mobility matrix, wedefine a linking matrix:L(Ω,α|{dot over (Λ)},β|c)=M(Ω,α|Λ,β|c)·{ P(Ω,α)• P(Λ,β)},  Eqn. 4.17

where 0≦L≦1.0.

If a template matching algorithm is used, then the product of maximizingcorrelations, via equation 4.12 and 4.13, defines the following linkingfactor:0≦U(Ω,α;Γ,β|c,Θ _(MAX))=Ĉ(Ω,α;c,Θ _(MAX))≦1.0,  Eqn. 4.18

which can be used in place of the product of ID vectors P(Ω,α)• P(Λ,β),i.e.,L(Ω,α|Γ,β|c)=M(Ω,α|Γ,β|c)·U(Ω,α;Λ,β|c,Θ_(MAX)).  Eqn. 4.19

In the limit that two frames are very close together, a good assumptionis that target pose is similar in both frames. If the pose relative tothe moving sensor platform is the same in both frames, one can usedirect correlation between the two raw signatures as defined in equation4.11 to derive a good similarity measure independent of class labels c.However, the pose relative to a moving sensor may not be invariant eventhough a target's orientation is constant in a fixed coordinate system.In the fixed coordinate system, the constant pose estimate is made bycomputing the average velocity vector between the two reports. Thesepose estimates are then converted into a specific orientations Θ_(spec.)in the sensor's coordinate system and then used to generate an HRRreference signatures Ψ_(Ref). (Θ_(Spec.), c) which are correlated withmeasured signatures via equation 4.12. The product of correspondingcorrelation measures from two frames define the following linkingfactor:V(Ω,α;Γ,β|c)=Ĉ(Ω,α;c,Θ_(α-Spec.))·Ĉ(Γ,β;c,Θ_(β-Spec.)).  Eqn. 4.20

The specific pose values, Θ_(Spec.), differentiate the product in 4.20from the product in equation 4.18 which utilizes poses, Θ_(MAX),associated with maximizing templates in the database.

In the limit that three frames are very close together, then one cancompute the average velocities between reports in adjacent frames anduse these velocities to project (forward and backwards) the position ofreports in the middle frame into the other two frames. The predictedextended locations are compared with the measured positions in the twoframes to derive a similarity measure. In general the measure would besome function of the separation between predicted {circumflex over (r(Ω;Γ.β) and measured r(Ω,α) locations; for example,D(Ω,α;Γ,β)=k·exp{({circumflex over ( r (Ω;Γ,β)−{right arrow over(r)}(Ω,α))²/σ²}.  Eqn. 4.21

This factor is not derived from ATR considerations; but uses a kinematicapproximation to generate relative linking likelihoods. Notice that theD matrix is similar to a predictive filtering technique, which isindependent of the class label c. This factor can be used to enhance thedefinition of the linking matrix or it can be use as a substitute factorfor the mobility matrix. In the case of high frame rates, it wouldsimplify computation to substitute the D matrix for the mobility matrix.

Several different components, U, V, D (equations 4-18, 4-20, 4.21) havebeen defined that could be used to redefine the linking matrix definedin equation 4.17. For example, for high frame rates, it would be moreconvenient to use V and D in the following definition:L(Ω,α|Γ,β|c)=D(Ω,α|Γ,β)·V(Ω,α;Γ,β|c,Θ_(α-Spec.),Θ_(β-Spec.)).

Also, if the pose relative to the sensor can be assumed to be the samein the two closely spaced frames, then the direct correlation factors,equation 4.11, could be the linking matrix components,L(Ω,α|Γ,β)=C(Ω,α|Γ,β),  Eqn. 4.23

In examples 4.22 and 4.23, the information that would be provided bymobility analysis is compensated for by higher frame rates. As frameseparations increase, it becomes necessary to replace these two factorswith an alternatives. Establishing the linking matrix between widelyseparated frames will depend heavily on the use of mobility analysis.However, it may be prudent to have several levels of mobility analysisto use for different frame separations. In this vein the D matrix wouldbe at the lowest level. When mobility analysis is unfeasible orresources limited, the linking matrix definition must rely solely on oneof the ATR factors discussed above.

It should be appreciated that different definitions of the linkingmatrix, dynamic models, and different ATR algorithms for generating theID vector may be used consistent with the present invention. Forinstance, different modes of computing the linking matrix areillustrated herein, which provides a good deal of flexibility informulating an HRR tracking problem.

4.4 Co-Moving Target Model

There is one final enhancement to the linking matrix that depends on thecontinuity of groups of co-moving targets. The concept is that a set oftargets moving together in one frame should tend to continue movingtogether in adjacent frames. We shall take a given report locatedoff-road at r(Ω,α) and define a neighborhood set of co-moving reports bythe following conditions:| r (Ω,α)− r (Ω,β)|≦δ,  Eqn. 4.24a|{dot over (ρ)}(Ω,α)−{dot over (ρ)}(Ω,β)|≦ε.  Eqn. 4.24b

The on-road conditions are 4.24a and the same speed tangent to the road:

$\begin{matrix}{{{{\frac{\overset{.}{r}\left( {\Omega,\alpha} \right)}{{{\overset{\_}{e}}_{\phi}\left( {\Omega,\alpha} \right)} \cdot {\overset{\_}{e}}_{ROAD\_\alpha}} - \frac{\overset{.}{r}\left( {\Omega,\beta} \right)}{{{\overset{\_}{e}}_{\phi}\left( {\Omega,\beta} \right)} \cdot {\overset{\_}{e}}_{ROAD\_\beta}}}} \leq ɛ},} & {{{Eqn}.\mspace{14mu} 4.24}c}\end{matrix}$

where the unit vector ē_(ROAD) _(—) _(X) is tangent to the road atreport X.

The analysis of co-moving targets is closely related to the problem ofdistinguishing multiple target reports in the same frame. Because thetime differences between reports in the same frame are relatively small,the constraints of closeness and matching speeds are similar toconstraints 4.24. Recall constraint 4.10, which assumed that the problemof multiple hits per frame has been eliminated elsewhere. This is ageneral problem for all tracking systems and a full discussion of itssolution is beyond the scoop of this paper. While leaving this subject,we point out that the correlation factor in equation 4.11 may be used totag two reports as belonging to different targets:C(Ω,α|Ω,β)≦threshold

different targets.

Let η(Ω,α|Ω,β)=1 if βε co-moving neighborhood of α=0 otherwise Eqn.4.25a

$\begin{matrix}{{\eta\left( {\Omega,\alpha} \right)} = {{\sum\limits_{\beta}{\eta\left( {\Omega,\left. \alpha \middle| \Omega \right.,\beta} \right)}} = {\alpha\text{-}{neighborhood}\mspace{14mu}{{size}.}}}} & {{{Eqn}.\mspace{14mu} 4.25}b}\end{matrix}$

Define the average ID for the α-neighborhood as

$\begin{matrix}\begin{matrix}{{{\hat{P}\left( {\left. C \middle| \Omega \right.,\alpha} \right)} = {{\sum\limits_{\beta,{\beta \neq \alpha}}{\frac{{\eta\left( {\Omega,\left. \alpha \middle| \Omega \right.,\beta} \right)} \cdot {P\left( {\left. C \middle| \Omega \right.,\beta} \right)}}{\eta\left( {\Omega,\alpha} \right)}{if}\mspace{14mu}{\eta\left( {\Omega,\alpha} \right)}}} > 0}},{C = 1},2,\ldots\mspace{14mu},N_{c}} \\{= {{0\mspace{14mu}{if}\mspace{14mu}{\eta\left( {\Omega,\alpha} \right)}} = 0}}\end{matrix} & {{Eqn}.\mspace{14mu} 4.26}\end{matrix}$

and the following linking functions:A(Ω,α|Γ,β)≡exp{−{η(Ω,α)−η(Γ,β)}²/σ²}≡exp{−Δη(Ω,α|Γ,β²)/σ²},  Eqn. 4.27

$\begin{matrix}{{{B\left( {\Omega,\left. \alpha \middle| \Gamma \right.,\beta} \right)} \equiv {\sum\limits_{c}{{\hat{P}\left( {\left. C \middle| \Omega \right.,\alpha} \right)} \cdot {\hat{P}\left( {\left. C \middle| \Gamma \right.,\beta} \right)} \cdot \left( {1 - {v \cdot {\delta(C)}}} \right)}} \equiv {{\overset{\underset{-}{\hat{}}}{P}\left( {\Omega,\alpha} \right)}\bullet{\overset{\underset{-}{\hat{}}}{P}\left( {\Gamma,\beta} \right)}}},} & {{Eqn}.\mspace{14mu} 4.28}\end{matrix}$N(Ω,α|Γ,β)={ω1+ω2·δ(η(Γ,β)}··A(Ω,α|Γ,β)+ω2·B(Ω,α|Γ,β),  Eqn. 4.29

where the positive weights satisfy ω1+ω2=1.0 and δ(X) is the deltafunction. Notice that A(Ω,α|Γ,β) is independent of the ATR derived IDvector and is a function of the mismatch between the different sizes ofco-moving neighborhoods. B(Ω,α|Γ,β) is a function of the ID vectors andequal to zero if the size of either of the two neighborhoods is equal tozero. The co-moving neighborhood function N(Ω,α|Γ,β) combines the twolinking functions and takes into account the case of null neighborhoods.

The co-moving neighborhood linking function is used to enhance thelinking matrix defined by equation 4.12. Let the enhanced definition begiven byL(Ω,α|Λ,β|C)≡W(Ω,α|Λ,β)·N(Ω,α|Λ,β)·M(Ω,α|Λ,β|C)  Eqn. 4.30a

or the alternative definition, which uses a weighted sumL(Ω,α|Λ,β|C)≡{v1·W(Ω,α|Λ,β)+v2·N(Ω,α|,Λ,β)}·M(Ω,α|Λ,β|C)  Eqn. 4.30b

where0≦L(M(Ω,α|Λ,β,C))≦1.0  Eqn. 4.31

and the positive weights satisfy, v1+v2=1.0.

4.5 Converting the Linking Matrix into a Compatibility Matrix

To covert a linking matrix into the basic compatibility matrix, we needto define a conditional probability distribution, P(C|j), over targetclasses (C) given the track index (j). Initially each unique track label(n) is anchored to a specific report (A,K) so it picks up the ID vectorof the anchoring report to define the required distributionP(c|n)≡P(c|κ,k).  Eqn. 4.32

It should also be noted that the new track (n) has an invariant labelweighting distribution in the anchoring frame ΛP(n|Λ,α)=δ _(ak),  Eqn. 4.33

where δ_(ak) is the Dirac delta function.

The conversion to the compatibility matrix is accomplished by usingP(c|j) to take a weighted average over the linking matrix,

$\begin{matrix}{{{R\left( {j,\Omega,\left. \alpha \middle| j \right.,\Lambda,\beta} \right)} \equiv {\sum\limits_{c = 0}^{N_{c}}{{{L\left( {\Omega,\left. \alpha \middle| \Lambda \right.,\left. \beta \middle| c \right.} \right)} \cdot {P\left( c \middle| j \right)}}\mspace{14mu}{for}\mspace{14mu}\alpha}}},{\beta > 0.}} & {{Eqn}.\mspace{14mu} 4.34}\end{matrix}$

The remaining components for the basic compatibility matrix involvemissing report labels (α or β=0). These components are the mostdifficult to define because constraints between missing reports andmeasured reports cannot be related to mobility analysis or ID vectors.An adaptive control definition for these components will be deferreduntil the other compatibility matrices are defined because it couplesall null labeled support functions.

4.6 Advanced Implementation Issues

Before concluding this section, it should be noted that for applicationswith large numbers of tracks, the compatibility matrix requires largeamounts of storage capacity. For a thousand targets, the basiccompatibility matrix has on the order of a billion components. In orderto reduce memory requirements, it is more efficient to store the linkingmatrix and use equation 4.34 to compute the coefficients as needed. Fora thousand track problem, the storage requirements would be severalorders of magnitude smaller.

The second implementation note also concerns equation 4.34 and is aspeculative suggestion for a new adaptive mechanism for futureexperimentation. After the initial phase of the relaxation process, asthe tracks become less ambiguous, each track may develop a set D ofdominant labels in each frame called anchors. One definition for adominant report for track j is given by:

If P(j,Ω,δ)>τ, then report (Ω,δ) is an anchor for track j. Eqn. 4.35

If τ>0.5 then there is only one dominant report for each track in eachframe; however, the same report can dominate two different tracks. Inorder to eliminate this ambiguity, we will design a negativecompatibility matrix in Section 5 to introduce competition betweentracks for the same report. Using the anchor set, one can define aweighted average ID vector over the set D:

$\begin{matrix}{\left\langle {P\left( c \middle| j \right)} \right\rangle = {\frac{\sum\limits_{\Omega,{\alpha \in D}}{{P\left( {j,\Omega,\alpha} \right)} \cdot {P\left( {\left. c \middle| \Omega \right.,\alpha} \right)}}}{\sum\limits_{\Omega,{\alpha \in D}}{P\left( {j,\Omega,\alpha} \right)}}.}} & {{Eqn}.\mspace{14mu} 4.36}\end{matrix}$

The average <P(c|j)> can be substituted for the anchor value ofP(c|j)=P(c|Ω₀,λ₀) in equation 4.34 that converts the linking matrix intothe compatibility matrix:

$\begin{matrix}{{{R\left( {j,\Omega,\left. \alpha \middle| j \right.,\Lambda,\beta} \right)} \equiv {\sum\limits_{c = 0}^{N_{c}}{{{L\left( {\Omega,\left. \alpha \middle| \Lambda \right.,\left. \beta \middle| c \right.} \right)} \cdot \left\langle {P\left( c \middle| j \right)} \right\rangle}\mspace{14mu}{for}\mspace{14mu}\alpha}}},{\beta > 0.}} & {{Eqn}.\mspace{14mu} 4.37}\end{matrix}$

This adaptive process introduces the concept of using track formation toimprove the ID vector which in turn modifies the compatibility matrixwhich then influences further track definition and so forth. For largeapplications in which one stores the linking matrix, this approach wouldrequire little additional computation.

5. Negative Compatibility Matrix (Lateral Inhibition)

The function of the negative compatibility matrix is to provide acompetitive pressure in the relaxation process that drives the labelingweight distribution to delta functions. In competitive neural networksthis process is called lateral inhibition or winner-take-all strategy.The idea is that if a certain report (Ω,α) has a large value of P(j,Ω,α)it will tend to suppress the weight associated with other reports havingsmaller components (P(i,Ω,α)<P(j,Ω,α)) so that P(j,Ω,α) increases at theexpense of P(i,Ω,α) for all i≠j. The negative compatibility matrix(denoted by R) is quite simple compared to the basic matrix discussed inthe previous section,

$\begin{matrix}{{\overset{\_}{R}\left( {i,\Omega,\left. \alpha \middle| j \right.,\Gamma,\beta} \right)} \equiv {{- \frac{\left( {1 - \delta_{ij}} \right)}{\left( {M_{\Omega} - 1} \right)}} \cdot \delta_{\Omega\Gamma} \cdot {\delta_{\alpha\beta}.}}} & {{Eqn}.\mspace{14mu} 5.1}\end{matrix}$

Note that R satisfies

${{- 1} \leq \frac{- 1}{\left( {M_{\Omega} - 1} \right)} \leq {\overset{\_}{R}\left( {i,\Omega,\left. \alpha \middle| j \right.,\Gamma,\beta} \right)} \leq 0},$

which is a necessary condition for theorems 3.1, 3.2, 3.3. This is trueeven if the normalizing parameter (M_(Ω)−1) were not included in thedefinition.

The corresponding support function is given by

$\begin{matrix}{{N\left( {i,\Omega,\alpha} \right)} \equiv {\sum\limits_{j,\Gamma,\beta}\;{{\overset{\_}{R}\left( {i,\Omega,\left. \alpha \middle| j \right.,\Gamma,\beta} \right)} \cdot {P\left( {j,\Gamma,\beta} \right)} \cdot {- {\sum\limits_{j = 1}^{M_{\Omega}}\;{\frac{\left( {1 - \delta_{ij}} \right)}{\left( {M_{\Omega} - 1} \right)} \cdot {{P\left( {j,\Omega,\alpha} \right)}.\mspace{20mu}{N\left( {i,\Omega,\alpha} \right)}}}}}}} \equiv {{- \frac{\left( {{\overset{\sim}{P}\left( {\Omega,\alpha} \right)} - {P\left( {i,\Omega,\alpha} \right)}} \right)}{\left( {M_{\Omega} - 1} \right)}}\mspace{14mu}{where}}} & {{{Eqn}.\mspace{14mu} 5.2}a} \\{\mspace{76mu}{{\sum\limits_{i = 1}^{M_{\Omega}}\;{N\left( {i,\Omega,\alpha} \right)}} = {{- {\sum\limits_{i = 1}^{M_{\Omega}}\;{P\left( {i,\Omega,\alpha} \right)}}} \equiv {{- {\overset{\sim}{P}\left( {\Omega,\alpha} \right)}}\mspace{14mu}{and}}}}} & {{Eqn}.\mspace{14mu} 5.3} \\{\mspace{70mu}{{\sum\limits_{\alpha = 1}^{N_{\Omega}}\;{N\left( {i,\Omega,\alpha} \right)}} = {- {1.0.}}}} & {{Eqn}.\mspace{14mu} 5.4}\end{matrix}$

For use in the next section, define

$\begin{matrix}{{{\hat{N}\left( {i,\Omega} \right)} \equiv {\sum\limits_{\alpha = 0}^{N_{\Omega}}\;{{P\left( {i,\Omega,\alpha} \right)} \cdot {N\left( {i,\Omega,\alpha} \right)}}}} = {- {\sum\limits_{\alpha = 0}^{N_{\Omega}}\;{{P\left( {i,\Omega,\alpha} \right)} \cdot {\frac{\left( {{\overset{\sim}{P}\left( {\Omega,\alpha} \right)} - {P\left( {i,\Omega,\alpha} \right)}} \right)}{\left( {M_{\Omega} - 1} \right)}.}}}}} & {{Eqn}.\mspace{14mu} 5.5}\end{matrix}$

All non-zero values of the support function are negative which tends toweaken the corresponding component of the assignment vector. Since thestronger components of the assignment vectors incur less negativesupport components they come out better in the competition based on thenegative compatibility matrix.

6. Combining Support Functions

Combining the three different support functions from above, we obtain aweighted sumS(i,Ω,α)=W ₁ ·H(i,Ω,α)+W ₂ ·V(i,Ω,α)+W ₃ N(i,Ω,α),  Eqn. 6.1

where

${\sum\limits_{Z = 1}^{3}\; W_{Z}} = 1.0$and 0≦W_(Z)≦1.0.

Since the basic and negative compatibility matrices satisfy theconditions of theorem 3.1 and are explicitly symmetric, we can apply thetheorems 3.1, 3.2, and 3.3 and use the gradient of the average localconsistency

$\begin{matrix}{{G\left( {i,\Omega,\alpha} \right)} = {\frac{\partial{A\left( \overset{\_}{P} \right)}}{\partial{P\left( {i,\Omega,\alpha} \right)}} = {{\sum\limits_{z = 1}^{3}\;{O_{Z} \cdot W_{Z} \cdot {S_{Z}\left( {i,\Omega,\alpha} \right)}}} = {{O_{1} \cdot W_{1} \cdot {H\left( {i,\Omega,\alpha} \right)}} + {O_{2} \cdot W_{2} \cdot {V\left( {i,\Omega,\alpha} \right)}} + {O_{3} \cdot W_{3} \cdot {N\left( {i,\Omega,\alpha} \right)}}}}}} & {{Eqn}.\mspace{14mu} 6.2}\end{matrix}$

where O₁=1 and O₂=O₃=2, to define the update rule:dP(i,Ω,α)=P(i,Ω,α)·dG(i,Ω,α)  Eqn. 6.3

wheredG(i,Ω,α)=G(i,Ω,α)−Ĝ(i,Ω),  Eqn. 6.4dG(i,Ω,α)=W ₁ ·dH(i,Ω,α)+2·W ₂ ·dV(i,Ω,α)+2·W ₃ ·dN(i,Ω,α),  Eqn. 6.5

and

$\begin{matrix}{{\hat{G}\left( {i,\Omega} \right)} = {{\sum\limits_{\alpha = 1}^{N_{\Omega}}\;{{P\left( {i,\Omega,\alpha} \right)} \cdot {G\left( {i,\Omega,\alpha} \right)}}} = {{W_{1} \cdot {\hat{H}\left( {i,\Omega} \right)}} + {2 \cdot W_{2} \cdot {\hat{V}\left( {i,\Omega} \right)}} + {2 \cdot W_{3} \cdot {{\hat{N}\left( {i,\Omega} \right)}.}}}}} & {{Eqn}.\mspace{14mu} 6.6}\end{matrix}$

The reason that O₁=1 is the definition of the local average consistencyA(P) spans the historic window in which the label assignment vectors aretreated as constants when computing the gradient.

We now present several useful theorems based on the update rule 6.3.

Theorem 6.1

Relationship P(j,QΩ,α)=δ_(αβ) is invariant under update rule 6.3.

Proof If P^(K)(j,Ω,α)=δ_(αβ),then

${\hat{G}\left( {j,\Omega} \right)} = {{\sum\limits_{\alpha = 1}^{N_{\Omega}}\;{{P^{K}\left( {j,\Omega,\alpha} \right)} \cdot {G^{K}\left( {j,\Omega,\alpha} \right)}}} = {G^{K}\left( {j,\Omega,\beta} \right)}}$dG ^(K)(j,Ω,α)≡(1−δ_(αβ))G ^(K)(j,Ω,α).

Substitution of dG into the update rule results in the followingidentitydP ^(K)(j,Ω,α)=P^(K)(j,Ω,α)·dG ^(K)(j,Ω,α)=δ_(αβ)·(1−δ_(αβ))·G^(K)(j,Ω,α)≡0 for all α.

Hence P^(K+1)(j,Ω,α)=P^(K)(j,Ω,α)=δ_(αβ). QED

Theorem 6.2

Relationship P(j,Ω,α)=0 for some α is invariant under update rule 6.3.

Proof: By definition,d ^(PK)(j,Ω,α)=P ^(K)(j,Ω,α)·dG ^(K)(j,Ω,α)=0 if P ^(K)(j,Ω,α)=0. HenceP ^(K+1)(j,Ω,α)=P ^(K)(j,Ω,α)=0.QED

In order to complete the formalism depending on the basic compatibilitymatrix, we now define the components of the basic matrix that havemissing report labels (zero labels). Let0≦R(i,Ω,0|i,γ,β)=R(i,Γ,β|i,Ω,0)=C _(Ω)≦1.0, for all Ω,Γ,i,β  Eqn. 6.7

so that active and historic support functions becomeV(i,Ω,0)=C _(Ω) or H(i,Ω,0)=C _(Ω),  Eqn. 6.8

where C_(Ω) is a variable parameter to be defined. The constraint0≦C_(Ω)≦1.0 is required by equation 4.2. There is ambiguity in zero—zerocomponents R(i,Ω,0|i,Γ,0) that could equal C_(Ω) or C_(Γ) by thedefinition 6.7. However, equations 6.8 are the more fundamentalassumptions in our model. Note, definition 6.7 is unambiguous forcomputing the support functions V(i,Ω,α) and H(i,Ω,α) when α>0.

The zero component of the gradient of the average local consistencyfunction (see definition 6.2) isG(i,Ω,0)=W ₁ ·H(i,Ω,0)+2·W ₂ ·V(i,Ω,0)+2·W ₃ ·N(i,Ω,0),G(i,Ω,0)=(W ₁+2·W ₂)·C_(Ω)+2·W ₃ ·N(i,Ω,0).  Eqn. 6.9

Summing over the track index i and noting equation 5.3 we obtain

$\begin{matrix}{{{\sum\limits_{i = 1}^{M_{\Omega}}\;{G\left( {i,\Omega,0} \right)}} \equiv {\overset{\sim}{G}\left( {\Omega,0} \right)}} = {{\left( {W_{1} + {2 \cdot W_{2}}} \right) \cdot M_{\Omega} \cdot C_{\Omega}} - {2 \cdot W_{3} \cdot {{\overset{\sim}{P}\left( {\Omega,0} \right)}.}}}} & {{Eqn}.\mspace{14mu} 6.10}\end{matrix}$

Intuitively, one can interpret {tilde over (P)}(Ω,0)/M_(Ω) as theprobability of a missing report in frame Ω. For the purpose of definingC_(Ω), an ideal value for this probability based on the differencebetween the number of active tracks and measured reports is estimated.Let the ideal value of the missing report probability be given by

$\begin{matrix}{= {\frac{1 + {\max\left( {{M_{\Omega} - N_{\Omega}},0} \right)}}{M_{\Omega}} > 0.}} & {{Eqn}.\mspace{14mu} 6.11}\end{matrix}$

Generally there should be more tracks than reports in a frame,M_(Ω)≧N_(Ω), otherwise one should generate new tracks. Let the value ofC_(Ω) be proportional to the ideal probability, which increases themissing report components of the compatibility matrix as the disparitybetween the number of tracks and reports increases:

$\begin{matrix}{{C_{\Omega} = {\frac{2 \cdot W_{3} \cdot}{\left( {W_{1} + {2 \cdot W_{2}}} \right)} > 0}},} & {{Eqn}.\mspace{14mu} 6.12}\end{matrix}$

where C_(Ω) must also be less than 1.0. Substituting C_(Ω) into equation6.10 yields

$\begin{matrix}{{{\frac{\overset{\sim}{G}\left( {\Omega,0} \right)}{M_{\Omega}} = {2 \cdot W_{3} \cdot \left( {- \frac{\overset{\sim}{P}\left( {\Omega,0} \right)}{M_{\Omega}}} \right)}},{where}}{{\overset{\sim}{P}\left( {\Omega,0} \right)} \equiv {\sum\limits_{i = 1}^{M_{\Omega}}\;{{P\left( {i,\Omega,0} \right)}\mspace{14mu}{and}}}}\mspace{14mu}{{\overset{\sim}{G}\left( {\Omega,0} \right)} \equiv {\sum\limits_{i = 1}^{M_{\Omega}}\;{{G\left( {i,\Omega,0} \right)}.}}}} & {{Eqn}.\mspace{14mu} 6.13}\end{matrix}$

Equation 6.13 relates the zero component averages of the gradient vectorand the assignment vector:G(Ω,0)/M _(Ω)>0 if {tilde over (P)}(Ω,0)/M _(Ω)<

,{tilde over (G)}(Ω,0)/M _(Ω)<0 if {tilde over (P)}(Ω,0)/M _(Ω)>

.  Eqn. 6.14

Equation 6.13 demonstrates that the definition of C_(Ω) is equivalent toa control mechanism that tries to balance the estimated probability ofmissing reports with the ideal probability dependent on the number oftracks and thereby influence the labeling distribution weights P(i,Ω,0)associated with missing reports. This control mechanism is a desirablefeature for modeling the difficult phenomenon of missing reports.

7. Initiating and Removing Tracks

The numbers of reports in each frame fluctuate for a variety of reasons.Those causing missing reports include the following: (1) targets stop orhave radial velocities relative to sensor which make them invisible, (2)targets leaving area of surveillance, (3) terrain occlusion(predictable) and other random masking of radar signals, and (4) sensorartifacts. Reasons that cause new target reports include: targetsstarting to move, entering the surveillance area, or emerging from amasking area. Also clutter or noise produces false target reports. Thesefluctuations present additional difficulties for track associationproblems such as the procedures for establishing estimates for missingreports discussed in the previous section. Means must be provided toadaptively adjust the number of tracks whose number cannot bepredetermined because of the above factors. Since the mathematicaltheory of RLA generally addresses problems in which the number ofobjects and labels are fixed, A-Track needed to explore new facets ofthe theory and provide a novel method for generating new tracks withinthe RLA model.

In the initial frame Θ, each sensor report is a unique target and isassigned to a new track label, M_(⊖)=N_(⊖) Each of the label weightdistributions is equal to a delta functionP(i,Θ,j)=δ_(ij),  Eqn. 7.1

where the number of track and report indices are equal in the initialframe. These initial assignments are invariant under the RLA updaterules as shown in theorem 6.1.

In frames other than the initial frame, A-Track uses the followingprocedure for generating and establishing a new track.

Criteria for Generating a New Track:

Certainly the condition, N_(Ω)>M_(Ω), of having more reports than tracksrequires one to generate new tracks as mentioned above. If one or morenew tracks have to be generated, one needs methods for anchoring eachnew track.

1. At any time, if

$\begin{matrix}{{{\overset{\sim}{P}\left( {\Omega,\lambda} \right)} \equiv {\sum\limits_{j = 1}^{M_{\Omega}}\;{P\left( {j,\Omega,\lambda} \right)}} < ɛ},} & {{Eqn}.\mspace{14mu} 7.2}\end{matrix}$

then a new track, with index I_(NEW), may be anchored to report (Ω,λ):P(I _(NEW),Ω,α)=δ_(αλ) and the ID vector P(C|I_(NEW))≡P(C|Ω,λ).

2. If a new track needs be generated and anchored to a report becauseN_(Ω)>M_(Ω), select that report corresponding to min_(α)({tilde over(P)}(Ω,α)).

Criteria 1 and 2 assert that the selected anchoring report (Ω,λ) isleast strongly associated with any track in the frame and therefore theappropriate candidate for establishing a new track.

After anchoring the new track to the candidate report (Ω,λ), one has thetask of establishing values for R(I_(NEW),Γ,β|I_(NEW),Ω,λ) andP(I_(NEW),Γ,α) for all Γ≠Ω. The required components of the compatibilitymatrix are computed as discussed in Section 4. No modification of theexisting components are necessary to accommodate new tracks. The ruleadopted in A-Track for initializing the distributions for the new trackin all active windows isR(I _(NEW),Γ,0|I_(NEW),Ω,λ)=C ₂,  Eqn. 7.3

$\begin{matrix}{{P\left( {I_{NEW},\Gamma,\beta} \right)} \equiv \frac{R\left( {I_{NEW},\Gamma,\left. \beta \middle| I_{NEW} \right.,\Omega,\lambda} \right)}{\sum\limits_{\beta = 0}^{N_{\Gamma}}\;{R\left( {I_{NEW},\Gamma,\left. \beta \middle| I_{NEW} \right.,\Omega,\lambda} \right)}} \geq 0} & {{Eqn}.\mspace{14mu} 7.4}\end{matrix}$

where C_(Ω) is defined in Section 6 (see Equation 6.12). This formulainitializes and normalizes the assignment vector for all reachablereports. Recall that any component initialized with a value of zero willalways retain this value. Thus all components of reachable reports aregiven a competitive chance with a weight proportional to the magnitudeof the new compatibility matrix.

Criteria for Track Removal:

A sequence of a fixed number of missing reports is a common criterionfor track removal in other tracking algorithms. The following is anothersimple criterion for dropping tracks from the active set of tracks:

$\begin{matrix}{C_{M} = {\frac{\sum\limits_{\Omega = 0}^{N_{F}}\;{P\left( {i,\Omega,0} \right)}}{N_{F} + 1} \geq \mu}} & {{Eqn}.\mspace{14mu} 7.5}\end{matrix}$

where 0≦C_(M)≦1.0 is a normalized missing report ratio and μ is athreshold value for track removal.

The number of non-ambiguous labels defined by

$\begin{matrix}{{{N(i)} = {{\sum\limits_{\Omega}{\sum\limits_{\alpha = 1}^{M_{\Omega}}{h\left( {{P\left( {i,\Omega,\alpha} \right)} - 1 + ɛ} \right)}}} \geq \tau}},} & {{Eqn}.\mspace{14mu} 7.6}\end{matrix}$

where h(x)=1 if x≧0 and =0 otherwise and ε is a small parameter, canalso be used to determine the activity of a track index. The criterionisIf N(i)≧τ

active track.  Eqn. 7.7

A sequence of missing sensor reports for a given track may beinterpreted in different ways. It could be a target that has stopped orexited the surveillance area or the track index could have beenassociated with a noisy report that is difficult to link with otherreports in adjacent frames. This distinction may have usefulintelligence value; however, it can be difficult to make. Suchdistinctions could be based on the pattern of missing reports. If,following the initialization of a track, one observes a series of stronglinks to observable reports in adjacent frames followed by a number ofmissing reports then one would assign a higher likelihood of the trackbeing a real target that has stopped or traveled away from thesurveillance area. Without an initial sequence of strong non-missingweights to verify the tracks authenticity, one tends to interpret thetrack as a noisy artifact. A large value of N(i) initially can be usedto establish a confidence measure for a track index being associatedwith a real target.

8. Third Order Compatibility Matrix

Thus far, the A-Track model does not exploit all the information one mayobtain from digital maps. For example, a digital terrain map allows oneto predict areas where one can expect to find missing reports. A digitalmap with environmental data would also allow one to forecast areas withhigher densities of false reports. This section illustrates how to use athird order compatibility matrix to exploit knowledge of predictableareas of sensor masking.

The A-Track masking model considers a triplet of report labelsassociated with the same track index in three separate frames. Thistriplet consists of a missing report temporally sandwiched between twoobservable reports. This model requires a third order compatibilitymatrix related to the masking support function by the following:

$\begin{matrix}{{M\left( {i,\Omega,\alpha} \right)} = {\sum\limits_{\Lambda,\Gamma}{D_{\Omega\Lambda\Gamma}{\sum\limits_{\beta,\gamma}{{R\left( {i,\Omega,{\alpha{{i,\Lambda,\beta}}i},\Gamma,\gamma} \right)} \cdot {P\left( {i,\Lambda,\gamma} \right)} \cdot {P\left( {i,\Gamma,\gamma} \right)}}}}}} & {{Eqn}.\mspace{14mu} 8.1} \\{where} & \; \\{{0 \leq {R\left( {i,\Omega,{\alpha{{i,\Lambda,\beta}}i},\Gamma,\gamma} \right)} \leq 1.0},} & {{Eqn}.\mspace{14mu} 8.2} \\{{D_{\Omega\Lambda\Gamma} \equiv {{\overset{\_}{\delta}}_{\Omega\Lambda} \cdot {\overset{\_}{\delta}}_{\Omega\Gamma} \cdot {\overset{\_}{\delta}}_{\Lambda\Gamma}} \equiv \frac{\left( {1 - \delta_{\Omega\Lambda}} \right) \cdot \left( {1 - \delta_{\Omega\Gamma}} \right) \cdot \left( {1 - \delta_{\Lambda\Gamma}} \right)}{F \cdot \left( {F - 1} \right)}},} & \; \\{{{\sum\limits_{\Lambda = 0}^{F}{\sum\limits_{\Gamma = 0}^{F}D_{\Omega\Lambda\Gamma}}} = 1.0},} & \;\end{matrix}$

δ _(XY)=1−δ_(XY) and δ_(XY) is the Dirac delta function. There should beno confusion between 0≦M (i,Ω,α)≦1.0 and the parameter M_(Ω) whichequals the number of tracks in frame Ω. The average local consistencyfor the compatibility matrix is

$\begin{matrix}{J_{M} = {\sum\limits_{i,\Omega,\alpha}{{M\left( {i,\Omega,\alpha} \right)} \cdot {P\left( {i,\Omega,\alpha} \right)}}}} & {{Eqn}.\mspace{14mu} 8.3} \\{or} & \; \\{J_{M} = {\sum\limits_{i,\Omega,\alpha}{\sum\limits_{j,\Lambda,\beta}{\sum\limits_{k,\Gamma,\gamma}{\delta_{ij} \cdot \delta_{ik} \cdot \delta_{jk} \cdot D_{\Omega\Lambda\Gamma} \cdot {R\left( {i,\Omega,{\alpha{{j,\Lambda,\beta}}k},\Gamma,\gamma} \right)} \cdot {P\left( {i,\Omega,\alpha} \right)} \cdot \left( {{P\left( {j,\Lambda,\beta} \right)} \cdot {{P\left( {k,\Gamma,\gamma} \right)}.}} \right.}}}}} & \;\end{matrix}$

The gradient is given by

$\begin{matrix}{\frac{\partial J_{M}}{\partial{P\left( {i,\Omega,\alpha} \right)}} = {\sum\limits_{\Lambda,\beta}{\sum\limits_{\Gamma,\gamma}{D_{\Omega\Lambda\Gamma} \cdot {P\left( {i,\Lambda,\beta} \right)} \cdot {P\left( {i,\Gamma,\gamma} \right)} \cdot {\left\{ {{R\left( {i,\Omega,{\alpha{{i,\Lambda,\beta}}i},\Gamma,\gamma} \right)} + {R\left( {i,\Lambda,{\beta{{i,\Omega,\alpha}}i},\Gamma,\gamma} \right)} + {R\left( {i,\Gamma,{\gamma{{i,\Lambda,\beta}}i},\Omega,\alpha} \right)}} \right\}.}}}}} & {{Eqn}.\mspace{14mu} 8.4}\end{matrix}$

We define the compatibility matrix by the following:R(i,Ω,α|i,Λ,β|i,Γ,γ)≡T(Ω,Λ,Γ))+T(Γ,Ω,Λ)+T(Λ,Γ,Ω)  Eqn. 8.5

whereT(ΩΛΓ)≡T(i,Ω,α|Λ°i,Γ,γ)≡ δ _(α0)·δ_(β0)· δ _(γ0)·A(Ω,α|Λ|Γ,γ)·B(Ω,Λ,Γ)·R(i,Ω,α|i,Λ,γ),T(ΓΩΛ)≡T(i,Γ,γ|Ωi,Λ,β)≡δ _(γ0)·δ_(α0) δ _(β0)·A(Γ,γ|Ω|Λ,β)·B(Γ,Ω,Λ)·R(i,Γ,γ|i,Λ,β),T(ΛΓΩ)≡T(i,Λ,β|Γ|i,Ω,α)≡ δ _(β0)·δ_(ε0)· δ _(α0)·A(Λ,β|Γ|Ω,α)·B(Λ,Γ,Ω)·R(i,Λ,β|i,Ω,α),B(Ω,Λ,Γ)≡H(Λ−Ω)·H(Γ−Λ), and H(x)=1 if x>0, =0 otherwise.

In the definition of function B, arguments of H are arithmeticdifferences of frame indices. The function B(Ω,Λ,Γ) imposes theconstraints that the three frame indices are different and that thecharacteristic time of each frame increases from left to right. Theproduct of three delta functions are non-zero only when the middle frameis associated with a missing report label and the other two frames withobservable reports (labels>0). The inclusion of the basic compatibilitymatrix term serves two functions. It is zero if the two observablereports are not connectable. Otherwise, it gives a relative likelihoodthat the observable reports are connected by track i.

The remaining factor, 0≦A(Ω,α|Λ|Γ.β)≦1.0, is related to the predictedmasking areas. To evaluate this function one needs to estimate thelocation of the missing report in the middle frame and determine whetheror not it lies within a predicted masking area. If the missing report'sestimated location is within the predicted masking area, then A isassigned a positive value and if not it is assigned a zero value (SeeFIG. 4). Hence, the T functions and ultimately the third ordercompatibility matrix will only be non-zero for those conditions thatsandwich a known masking area between two connectable observations.

Mobility analysis may be used to estimate a missing report position fromthe two observable reports conditioned on them being connectable. Thesimplest estimate for the position of the missing report is to assume astraight-line trajectory between the observable reports:

$\begin{matrix}{{\overset{\_}{r}(\Lambda)} \approx {{\overset{\_}{r}\left( {\Omega,\alpha} \right)} + {\frac{{\overset{\sim}{t}(\Lambda)} - {t\left( {\Omega,\alpha} \right)}}{{t\left( {\Gamma,\beta} \right)} - {t\left( {\Omega,\alpha} \right)}} \cdot \left\{ {{\overset{\_}{r}\left( {\Gamma,\gamma} \right)} - {\overset{\_}{r}\left( {\Omega,\alpha} \right)}} \right\}}}} & {{Eqn}.\mspace{14mu} 8.7}\end{matrix}$

where {tilde over (τ)}(Λ) is a time characteristic of the Λ frame. Ifthe missing report region is large, then even this rough estimate may begood enough to establish non-zero constraint values forA(Ω,α|Λ|Γ.β);i.e.,If r (Λ)ε Masked Region, then A=1;A=0 otherwise.  Eqn. 8.8

This simplified algorithm can be replaced to achieve better estimates.For example, the straight-line trajectory assumption in 8.7 needs to bechanged to accommodate on-road trajectories for curving roads. Also,mobility analysis can determine multiple possible trajectories betweenthe two anchoring reports with some fraction going through the predictedmasking area in the mid-frame and others missing it. Possibletrajectories missing the masking areas are identified by the followingconstraint on the basic compatibility matricesR(i,Ω,α|i,Λ,β)·R(i,Λ,β|i,Γ,γ)>0,  Eqn. 8.9

where (Λ,β) is a linking observable report in the mid-frame. Themagnitude of A(Ω,α|Λ|Γ.β) is inversely proportional to the number oflinking reports satisfying this constraint 8.8 (see FIG. 5). Elaborationof an algorithm for computing A is beyond the scope of this discussion;however, it is easy to appreciate that it would be computed as aby-product of the module that does the mobility analysis.

The gradient of the average local consistency is given by

$\begin{matrix}{\frac{\partial J_{M}}{\partial{P\left( {i,\Omega,\alpha} \right)}} = {\sum\limits_{\Lambda,\beta}{\sum\limits_{\Gamma,\gamma}{{D({\Omega\Lambda\Gamma})} \cdot {P\left( {i,\Lambda,\beta} \right)} \cdot {P\left( {i,\Gamma,\gamma} \right)} \cdot \left\{ {{T({\Omega\Lambda\Gamma})} + {T({\Gamma\Omega\Lambda})} + {T({\Lambda\Gamma\Omega})} + {T({\Lambda\Omega\Gamma})} + {T({\Gamma\Lambda\Omega})} + {T({\Omega\Gamma\Lambda})} + {T({\Gamma\Lambda\Omega})} + {T({\Omega\Gamma\Lambda})} + {T({\Lambda\Omega\Gamma})}} \right\}}}}} & {{Eqn}.\mspace{14mu} 8.9} \\{or} & \; \\{\frac{\partial J_{M}}{\partial{P\left( {i,\Omega,\alpha} \right)}} = {\sum\limits_{\Lambda,\beta}{\sum\limits_{\Gamma,\gamma}{{D({\Omega\Lambda\Gamma})} \cdot {P\left( {i,\Lambda,\beta} \right)} \cdot {P\left( {i,\Gamma,\gamma} \right)} \cdot {\left\{ {{T({\Omega\Lambda\Gamma})} + {2 \cdot {T({\Omega\Gamma\Lambda})}} + {T({\Gamma\Omega\Lambda})} + {2 \cdot {T({\Lambda\Omega\Gamma})}} + {T({\Lambda\Gamma\Omega})} + {2 \cdot {T({\Gamma\Lambda\Omega})}}} \right\}.}}}}} & \;\end{matrix}$

Because of the symmetric product P(i,Λ,β)·P(i,Γ,γ) in the summation, onecan interchange the dummy indices sets (Λ,β) and (Γ,γ) so that

$\begin{matrix}{\frac{\partial J_{M}}{\partial{P\left( {i,\Omega,\alpha} \right)}} = {3 \cdot {\sum\limits_{\Lambda,\beta}{\sum\limits_{\Gamma,\gamma}{{D({\Omega\Lambda\Gamma})} \cdot {P\left( {i,\Lambda,\beta} \right)} \cdot {P\left( {1,\Gamma,\gamma} \right)} \cdot {\left\{ {{T({\Omega\Lambda\Gamma})} + {T({\Gamma\Omega\Lambda})} + {T({\Lambda\Gamma\Omega})}} \right\}.}}}}}} & {{Eqn}.\mspace{14mu} 8.10}\end{matrix}$

Substituting definition 8.5 into equation 8.10, shows that the gradientis three times the support function:

$\begin{matrix}\begin{matrix}{\frac{\partial J_{M}}{\partial{P\left( {i,\Omega,\alpha} \right)}} = {3 \cdot {\sum\limits_{\Lambda,\beta}{\sum\limits_{\Gamma,\gamma}{{D({\Omega\Lambda\Gamma})} \cdot {P\left( {i,\Lambda,\beta} \right)} \cdot {P\left( {i,\Gamma,\gamma} \right)} \cdot}}}}} \\{R\left( {i,\Omega,{\alpha{{i,\Lambda,\beta}}i},\Gamma,\gamma} \right)} \\{= {3 \cdot {{M\left( {i,\Omega,\alpha} \right)}.}}}\end{matrix} & {{Eqn}.\mspace{14mu} 8.11}\end{matrix}$

The compatibility matrix as defined by equation 8.5 does not satisfy the‘explicit’ symmetry constraint defined by equation 2.14; however, it has‘implicit’ symmetry which also yields the format of equation 8.11. Thedefinition of the 3^(rd) order compatibility matrix can be modified tohave explicit symmetry but this is not necessary since the implicitsymmetry is sufficient and more efficient. Recall that having thesymmetry, explicitly or implicitly, simplifies the computation of thegradient of the average local consistency used in the RLA updatingequation 3.23.

The equations for combining support functions given in Section 6 need tobe modified to include an additional function M (i,Ω,α). Some of themodifications necessary to define the total gradient are given in thefollowing summary:

$\begin{matrix}\begin{matrix}{{{G\left( {i,\Omega,\alpha} \right)} = {\frac{\partial{A\left( \overset{\_}{P} \right)}}{\partial{P\left( {i,\Omega,\alpha} \right)}} = {\sum\limits_{z = 1}^{4}{O_{Z} \cdot W_{Z} \cdot {S_{Z}\left( {i,\Omega,\alpha} \right)}}}}},} \\{= {{W_{1} \cdot {H\left( {i,\Omega,\alpha} \right)}} + {2 \cdot W_{2} \cdot {V\left( {i,\Omega,\alpha} \right)}} +}} \\{{{2 \cdot W_{3} \cdot {N\left( {i,\Omega,\alpha} \right)}} + {3 \cdot W_{4} \cdot {M\left( {i,\Omega,\alpha} \right)}}},}\end{matrix} & {{{Eqn}.\mspace{14mu} 6.2}a}\end{matrix}$dP(i,Ω,α)=P(i,Ω,α)·dG(i,Ω,α),  Eqn. 6.3dG(i,Ω,α)=G(i,Ω,α)−Ĝ(i,Ω),  Eqn. 6.4dG(i,Ω)=W ₁ ·dH(i,Ω,α)+2·W ₂ ·dV(i,Ω,α)+2·W ₃ ·dN(i,Ω,α)+3·W ₄·dM(i,Ω,α),  Eqn. 6.5a

$\begin{matrix}{{\hat{G}\left( {i,\Omega} \right)} = {{\sum\limits_{\alpha = 1}^{N_{\Omega}}{{P\left( {i,\Omega,\alpha} \right)} \cdot {G\left( {i,\Omega,\alpha} \right)}}} = {{W_{1} \cdot {\hat{H}\left( {i,\Omega} \right)}} + {2 \cdot W_{2} \cdot {\hat{V}\left( {i,\Omega} \right)}} + {2 \cdot W_{3} \cdot {\hat{N}\left( {i,\Omega} \right)}} + {3 \cdot W_{4} \cdot {{\hat{M}\left( {i,\Omega} \right)}.}}}}} & {{{Eqn}.\mspace{14mu} 6.6}a}\end{matrix}$

The modification to equations 6.13 is

$\begin{matrix}{{\frac{\overset{\sim}{G}\left( {\Omega,0} \right)}{M_{\Omega}} = {{W_{3} \cdot \left( {- \frac{\overset{\sim}{P}\left( {\Omega,0} \right)}{M_{\Omega}}} \right)} + {3 \cdot W_{4} \cdot \frac{\overset{\sim}{M}\left( {\Omega,0} \right)}{M_{\Omega}}}}},} & {{Eqn}.\mspace{14mu} 6.13}\end{matrix}$

where the additional support function satisfies

$\begin{matrix}{0 \leq {M\left( {i,\Omega,0} \right)} \leq {1.0\mspace{14mu}{and}\mspace{14mu} 0} \leq \frac{\overset{\sim}{M}\left( {\Omega,0} \right)}{M_{\Omega}} \equiv \frac{\sum\limits_{i = 1}^{M_{\Omega}}{M\left( {i,\Omega,0} \right)}}{M_{\Omega}} \leq {1.0.}} & {{Eqn}.\mspace{14mu} 8.12}\end{matrix}$

9. Operation of A-Track

By virtue of the foregoing, A-Track represents a new method for thedevelopment of tracking algorithms and is characterized by the use offormalism from relaxation labeling algorithms and flexibility formodeling diversified operational conditions. It may be driven by ATRalgorithms and mobility analysis. The HRR problem used to illustratethis flexibility demonstrates how different modalities can accommodatedata collection with wide variation in the temporal separation offrames. For large separation, the number of possible links increases andplaces greater demands on the association problem. This problem iscomplicated further by sensor artifacts, missing reports, and targetsleaving and entering the surveillance area. In order to make correctassociations as entropy increases, more of the available knowledge needsto be exploited.

The use of compatibility matrices as defined in relaxation labelingtheory is shown to be a convenient and easily adaptable formalism tomodel and integrate ATR inputs with mobility analysis. In order tofacilitate modeling of the HRR application it has been necessary toenhance the standard RLA outlined in Section 2. A number of new theoremsare presented in Section 3 that serve as the bases for severalimprovements. First, new theorems are presented which demonstrate thatby limiting the dynamic range of the compatibility matrix withoutsacrificing its symmetry, one can utilize the Chen-Luh simplified updaterule. Maintaining the symmetric properties simplifies computation of thegradient of the average local consistency, which is shown to be greaterthan or equal to zero, thus insuring local convergence of the model'sobjective function. Also, these new algorithms facilitate theintegration of four different support functions derived from threedifferent compatibility matrices defined in the A-Track HRR model. Thisset of symmetric matrices has different dynamic ranges and generatesdifferent support functions with different orders. In the resultingformalism, the gradient of the objective function equals a weightedaverage of the four support functions.

Ideally invariant features should be used as ATR inputs into the linkingmatrices computation; however, in general, no invariant HRR features areknown, other than the target ID vector. Recall that HRR signatures arevery sensitive to changes in target and sensor orientation. The IDvector, which is a generic output of ATR algorithms, can be easilyintegrated into the computation of the linking matrices. HRR signaturetemplates were discussed from the standpoint of generating signaturecorrelation measures and pose estimation. Not included in this A-Trackmodel are sensors other than HRR, such as ELINT, which provide poseinvariant features that correlate over large time intervals.

Mobility analysis is mainly concerned with using physical constraints toestablish relative likelihoods for possible connections between pairs ofobserved reports. The basic compatibility matrix for the HRR modelpresented in Section 4 demonstrates how to combine mobility analysiswith ATR input for widely separated frames and how to use predictivefiltering between closely separated frames. Coefficients in this secondorder matrix link two observed sensor reports. Conceptually moredifficult to define are the coefficients associated with missingreports. A solution is given in Section 6 in which all the missingreport coefficients in a frame are assigned the same parametric value.These parameters act as control variables that are adjusted eachiterative cycle in a manner that effects a balance between the number oftracks and number of reports in each frame. The basic matrix defines asecond order active support function linking only frames in the activesliding window and a first order historic support function that linksone active frame and recent frames in an adjacent inactive window.

In order to promote competition between report labels for each of thetracks, a second compatibility matrix with negative coefficients waspresented in Section 5. This competition tends to eliminate theambiguous case where several different tracks are all assigned largeassignment weights corresponding to the same report. Recall that the sumof assignment weights for the same report across all tracks is notnormalized as is the sum over alpha index in the assignment vectorP(i,Ω,α). Derived from basic and negative compatibility matrices arethree support functions whose unification is discussed in Section 6.

Algorithms for deleting old tracks and for generating and initializingnew tracks are given in Section 7. In most RLA applications the numberof objects and labels are fixed parameters. The number of tracks neededis a variable of the optimization process and varies as new frames areadded. Also the number of report labels varies from frame to frame. Inthis formalism the addition of a new track does not effect previouscomputations. With each new track the corresponding linking matrices arecomputed between the new track anchored to a specific report and allother frames. This in turn is used to initialize the assignment weightsfor the new track in all other frames in the active window. Alsodiscussed are methods for determining which disappearing tracks are realand which are noisy artifacts.

An additional compatibility matrix is presented in Section 8 in order toexploit knowledge of predicted masking areas not included in the basicmatrix. This enhancement demonstrates that capturing additionalconstraints may require a new compatibility matrix rather thanmodification of factors in an existing matrix. This new case required athird order compatibility matrix to span three different frames. Indesigning the masking model, the third order matrix utilized the secondorder basic matrix as a multiplicative factor in its definition. Thethird order matrix has implicit symmetry rather than the explicitsymmetry generally discussed in the literature. It is shown thatimplicit symmetry is efficient while allowing the gradient to beproportional to the support function. At the end of Section 8, themethodology for combining the fourth support function with the threeprevious support functions is given.

10. Implementation

FIG. 5 depicts an illustrative sequence of operations or procedure 100of an ATR assisted tracking algorithm implemented by the A-Track system10 described above. The formulations presented above are used to developsymmetric compatibility matrices so that equations for updatingassignment vectors are greatly simplified.

One or more sensors provide reports (block 102) that are grouped into anew frame as the moving windows are updated (block 104). The processingalgorithm may advantageously adapt to the temporal frame separation(block 106). For instance, certain targets may lend themselves totraditional tracking techniques with small temporal frame separations.Thus, the basic compatibility matrices are modified for the appropriatetemporal separation, an example of which is described below with regardto FIG. 6.

Thereafter, the reports for the new frame are modeled and analyzed toproduct a linking matrix (block 108). The target represented by thereport is linked with previous reports assigned to a track. The modelingtakes advantage of known constraints, such as from a target database 110and a road/terrain database 112. Thereby, certain links between reportsmay be precluded, simplifying the linking.

The linking matrix is then organized and assembled into compatibilitymatrices (block 114). Then, a relaxation labeling process is performedto update an assignment vector (block 116). Specifically, thecompatibility matrices are converted into support functions that arethen used to modulate the assignment vector P^(K). The relaxationlabeling loop is an iterative process for optimizing the objectivefunction the local average consistency function A (see equation 3.32).

The results of the assignments may indicate the need for adding ordeleting a track in the active window (block 118). Adding a new trackmay be warranted when an unassigned report remains after all tracks areassigned. A track may be deleted when a criteria is met, such as acertain number of frames are processed without assigning a report to thetrack or the modeling indicates that the track has departed thesurveillance area.

An optional coupling mechanism improving the conditional probabilityvector by taking a weighted average of the ID vectors (“target vectors”)with the assignment vectors (see equation 4.36). The conditionalprobability vector is improved by communicating the updated assignmentvector (block 120) back to block 108 for the next iteration of theprocess. During which, ID vectors in different frames are linked (block122). Without this option, the initial conditional probability vectorwould not change. The optional method for computing a weighted averageconditional probability vector, described in Section 4.6 (AdvancedImplementation Issues), provides a new approach for coupling the targetID vector derived from ATR analysis with the tracking information in theform of the assignment vector derived from the relaxation labelingalgorithm. A weighted average of the ID vector for a fixed trajectoryindex is taken across multiple frames using the assignment vectors. Thiscoupling is exciting because it allows the ATR and tracking componentsto influence and improve one another in the iterative relaxation loop.This is a unique feature of A-Track system that uses improved trackdefinition via assignment vectors to improve the conditional probabilityvector (see equation 4.36) and then uses it to modulate theCompatibility Matrix (see equation 4.37) so as to further improve theassignment vectors. This interaction between ATR and tracking continuesin an interactive process. The marriage of ATR and Tracking algorithmswas one of the main motivations in the derivation of A-TRACK.

Procedure 100 concludes each iteration with a trajectory descriptionoutput (block 122) that may be used by user for tracking vehicles and isalso fed back for use in the next iteration of the procedure.

FIGS. 6 and 7 depicts a sequence of operations or procedure 200 thatdescribes a more detailed implementation of an ATR assisted trackingsystem. In particular, procedure 200 illustrates adapting the basiccompatibility matrices to temporal spacing, using a compatibility matrixfor masking, using a negative compatibility matrix to avoid ambiguitybetween reports, and coupling target ID vectors with the assignmentvector. In addition, the above-described formulations arecross-referenced to the procedure 200.

Procedure 200 as depicted in FIG. 6 generally describes a process ofcomputing the basic compatibility matrices (blocks 202–226), which tendsto be complicated with increased frame temporal spacing.

Each iteration of the procedure 200 begins with receipt of frame reportsfrom one or more sensors (block 202). These sensor inputs include targetstate vector and target signatures. The moving windows (active andhistoric) are updated with the frame reports (block 204). Thereafter adetermination is made as to whether the temporal frame separation issmall (block 206). If small, then ATR target signature correlation isperformed (see Sect. 4.3 above) (block 208), referencing predeterminedinformation in a terrain/target characteristics database 210. Thisprepared database stores digital terrain maps used to estimate themasking areas and has access to navigational data as required to locatepositions on the digital map.

The results of which produce ATR output terms for a linking matrix (seeSect. 4) (block 212). In particular, a conditional probability vectorP(c|j), where c and j are the class and track indices respectively, areinitially defined to be equal to the ID target vector of a specificreport to which the new track is anchored in a specific frame (see Eqn.4.32). The assignment vector for a new track has an invariant labelweighting distribution in the anchoring frame (see Eqn. 4.33). Thisconditional probability vector is used to convert the linking matrixinto basic compatibility matrices for the active and historic windows(see Eqn. 4.34), as will be described below with respect to block 226.Predictive matching technique and target signature correlation (seeSect. 4.3) is performed on the ID target vector are then performed(block 214).

If the temporal frame separation was not small in block 206, then alinking matrix is produced taking advantage of information that canstill be gleaned through ATR analysis (see Sect. 4.3) (block 216). Inparticular, ATR analysis may advantageously provide a number of factors(generally multiplicative factors) to enhance the definition of thelinking matrix. For instance, the dot product of ID vectors (see Eqns.4.14–4.17 and related discussion); template matching factors (see Eqns.4,18–4.19) that also provide some target pose information; andcorrelation measures between target signatures in two frames (see Eqns.4.20, 4.22, and 4.23). These ATR analyses advantageously referenceterrain/target characteristics database 210 is assist in targetidentification that can constrain or determine a likelihood of variousassignments. Thereby ATR output terms for the linking matrix areproduced (see Sect. 4) (block 218), as described above for block 212.

Mobility analysis (see Sect. 4.2) and co-moving target model (see Sect.4.4) are thereafter performed, thereby advantageously refining theprobabilities. Mobility analysis provides a number of factors (generallymultiplicative) for defining the linking matrix and factors forpredicting masking areas (see Eqns. 8.7, 8.8 and Sect. 8 in general).

With reference to block 224, an optional coupling mechanism improves theconditional probability vector by taking a weighted average of the IDvectors (“target vectors”) (See Eqn. 4.32) with the assignment vectors(see Eqn. 4.36). The conditional probability vector is improved bycommunicating the updated assignment vector from the previous iterationof the process. During which, ID vectors in different frames are linked.Without this option, the initial conditional probability vector wouldnot change. The optional method for computing a weighted averageconditional probability vector, described in Section 4.6 (AdvancedImplementation Issues), provides a new approach for coupling the targetID vector derived from ATR analysis with the tracking information in theform of the assignment vector derived from the relaxation labelingalgorithm. A weighted average of the ID vector for a fixed trajectoryindex is taken across multiple frames using the assignment vectors. Thiscoupling is exciting because it allows the ATR and tracking componentsto influence and improve one another in the iterative relaxation loop.This is a unique feature of A-Track system that uses improved trackdefinition via assignment vectors to improve the conditional probabilityvector (see equation 4.36) and then uses it to modulate theCompatibility Matrix (see equation 4.37) so as to further improve theassignment vectors. This interaction between ATR and tracking continuesin an interactive process. The marriage of ATR and Tracking algorithmswas one of the main motivations in the derivation of the A-Track system.

After either block 214 or block 218, then the appropriate linking matrixL is adaptively assembled (see Sect. 4) (block 214) with reference tothe terrain/target characteristics database 210. The linking matrix L isan efficient way to store the assembled information necessary to computethe compatibility matrices. Then the basic compatibility matrices(R_(V), R_(H)) for the active and historic windows are formed by takinga weighted average of the linking matrix using the ID vector as shown inEqns. 4.32, 4.33.

In addition to the basic compatibility matrices R_(V), R_(H), FIG. 6also depicts enhancing the analysis with a predicting masking areas. Inparticular, the frame reports from block 202 and the terrain/targetcharacteristics database 210 are used in block 228 to estimate maskingareas. Then, a third-order compatibility matrix R_(M) (see Eqn. 8.5) isformed (block 230).

As a further addition to the basic compatibility matrices R_(V), R_(H),FIG. 6 also depicts forming a negative compatibility matrix (see Eqn.5.1) that adds competition between assignment vectors (block 232) toassist in the convergence to a single trajectory description during therelaxation label process.

With particular reference to FIG. 7, the remainder of the procedure 200is depicted, and in particular the relaxation label algorithm (RLA)process of blocks 238–246. This RLA process begins computing a supportfunction for each compatibility matrix that is the weighted average ofthe respective compatibility matrix with the assignment vector. Thesesupport functions are iteratively optimized in the RLA loop. As thesupport increases for associating a given track with a given targetreport, the assignment tends to follow this increase in the supportfunction. In a similar manner, the assignment vector will decrease forassociations that have smaller value of the support function. Inparticular, the masking support function M is computed (see Eqn. 8.1)(block 238). The basic support function V for the active window (seeEqn. 4.5) is computed (block 234). The basic support function H for thehistoric window is computed (see Eqn. 4.8) (block 236). The negativesupport function N is computed (see Eqn. 5.2) (block 240).

The objective function to be optimized is formed by averaging the localconsistency defined initially in Eqn. 2.11 and re-defined in Eqn. 3.10.(see Eqn. 6.1). If the RLA process is used to optimize the objectivefunction, then a simple formula may be used in terms of the gradient andaverage gradient of the objective function (see Eqn. 6.1–6.5, 6.6, 6.6A,6.13) (block 244).

The assignments vectors are updated in block 246 (see Eqns. 6.1, 6.2,6.3–6.6) with the updated assignment vector fed back to block 242 forrecalculating the support functions. The RLA loop of blocks 242–246 mayrepeat until ideally the assignment vectors relax (converge) to deltafunctions with a single component equal to 1 and the rest 0 (i.e., theassignment vectors become trajectory description outputs).

When the RLA loop is complete for this iteration of the overallprocedure 200, tracks are added or deleted as appropriate (see Sect. 7)(block 248) and then the trajectory description output is provided(block 250).

While the present invention has been illustrated by description ofseveral embodiments and while the illustrative embodiments have beendescribed in considerable detail, it is not the intention of theapplicant to restrict or in any way limit the scope of the appendedclaims to such detail. Additional advantages and modifications mayreadily appear to those skilled in the art. For example, A-Track is notlimited to the HRR application used in this paper to illustrate it. Italso provides the basis for modeling other tracking problems utilizingdifferent types of sensors and image analysis problems that need tointegrate a complicated suite of constraints. In ATR-driven trackingproblems, the full interaction of the dynamics between forming tracksand improving target ID vectors may include continuing to modify thebasic compatibility matrix by using a weighted average of the ID vectorsis outlined at the end of Section 4.

TABLE 1 Fundamental differences between two tracking applications usingdifferent relaxation algorithms. Lagrangian Relaxation RelaxationLabeling Algorithm Algorithm (Probabilistic Labeling) Variables Set ofcomplete tracks Set of weighted labeling assignments Z_(ijk...) ∈ {0, 1}binary variables P(i, Ω, α) real functions satisfying satisfyingconstraint equations. probability axioms. Prior Set of cost values foreach Set of flexible compatibility matrices Knowledge trajectory toincorporate all for incorporating various knowledge knowledge. aspects.Link all frames in window. Link small set of frames in window. ObjectiveTotal cost to be minimize Average local consistency to be FunctionLinear function of Z_(ijk...) maximized Non-linear function of P(i, Ω,α) Relaxation Seeks to satisfy constraint Modifies P(i, Ω, α) tomaximize Process equations while minimizing objective function definedby total cost of all trajectories compatibility matrices. utilized toconnect sensor reports. Constraints Set of complicating constraint Noconstraint equations required. Equations equations to insure precisetrajectories Output Complete and unique set of un- Set of weightedlabeling assignments ambiguous tracks specified by P(i, Ω, α) satisfyingprobability axioms Z_(ijk...) which are ambiguous track descriptorsPrecision Precise tracks Fuzzy track descriptions. Non-ambiguousdefinitions of May result in strictly-ambiguous full set completetracks. definitions of track descriptors P(i, Ω, α). Tolerates No YesAmbiguity Incurs complexity to avert Incurs ambiguity to avertcomplexity. ambiguity. Object. Set of constraint equations Set ofweighted labeling assignments Relaxed P(i, Ω, α) Complexity HIGH, due tocomplicate LOW, with simple set of updating multiple recasting ofproblem equations formulation. p^(K+l)(i, Ω, α) ≡ P^(K)(i, Ω, α) +dP^(K)(i, Ω, α) Special ATR None Factorization of compatibility matricesCoupling into linking matrices and ATR Mechanisms dependent conditionalprobability distribution. Natural Not evident. Averaged conditionalprobability feedback to distribution ATR systemfor dynamiccoupling oftwosystems$\left\langle {P\left( {c❘j} \right)} \right\rangle = \frac{\sum\limits_{\Omega,\alpha}{{P\left( {j,\Omega,\alpha} \right)} \cdot {P\left( {{c❘\Omega},\alpha} \right)}}}{\sum\limits_{\Omega,\alpha}{P\left( {j,\Omega,\alpha} \right)}}$is useful for upgrading both ATRACT and ATR systems and combinesvariables originating in both systems. Recasting of Required to recastbinary No recasting required. problem partition problem into differentformulation forms in order to optimize. Disambiguation Not required ifconstraint Generally required, but resulting Of relaxation equations aresatisfied. ambiguous assignments may be results changed into preciseassignments.

As described in Table 1, above, the generally-known LagrangianRelaxation Algorithms have several characteristics. First, Basicassignment variable: A binary assignment variable Z_(ijk) . . . ε{0,1}defines a specific trajectory with each index specifying a unique sensorreport or missing report for each frame in the moving window. In otherwords, the existence of a hypothetical possible trajectory consisting ofone assignment per frame for every frame is represented by the binaryvariable Z_(ijk) . . . .

Second, Prior knowledge: All prior knowledge is represented as a scoreor cost function which associates a cost value C_(ijk) . . . for everyprecisely defined trajectory designated by a binary assignment variableZ_(ijk) . . . .

Third, Objective function: Basic objective function is total cost to beminimized,

${Cost} = {\sum\limits_{i}{\sum\limits_{j}{\sum\limits_{\ldots}{C_{ijk\ldots} \cdot Z_{ijk\ldots}}}}}$However the binary assignments must satisfy a system of constraints thatrequire that all reports in a given frame be associated with one trackonly and that each track in a given frame be associated with only onesensor report or missing sensor report. The simultaneous satisfaction ofthis system of constraints complicates the process of minimizing thebasic objective function. This binary partitioning problem could besolved by enumerating all possible assignments to find the optimalassignment satisfying all constraints. However, this rapidly becomes animpractical approach for larger problems because of the vast number ofpossibilities. Therefore, other optimization techniques are employedsuch as a Lagrangian relaxation technique which needs to recast thecomplicated binary partitioning problem into a continuous portioningproblem from which a binary solution can ultimately be recovered.

It is the constraint equations that are relaxed by the Lagrangianrelaxation algorithm.

Fourth, Outputs: Set of non-zero binary trajectory assignments Z_(ijk) .. . ε{0,1} defining a set of precise trajectories, i.e. each reportlabel assigned to one track and each track assigned to one report labelor the missing report label.

By contrast, the inventive Relaxation Labeling Algorithm has a number ofcharacteristics with a number of advantages. First, Basic assignmentvariable: The basic assignment variables are called the weight labelassignments P(i,Ω,α) which satisfies the probability axioms,

${0 \leq {P\left( {i,\Omega,\alpha} \right)} \leq {1\mspace{14mu}{and}\mspace{14mu}{\sum\limits_{\alpha}{P\left( {i,\Omega,\alpha} \right)}}}},$and specify the relative likelihood that sensor report labels (Ω,α)being associated with the track object designated by (i) in frame (Ω).The system of weighted labeling assignments is generally astrictly-ambiguous labeling which does not specify a system of precisetrajectories.

Second, Prior knowledge: Different aspects of prior knowledge arerepresented by a system of varied-ordered compatibility functionslinking a finite number of reports in various frames. Differentcompatibility functions link a report to other reports in the movingactive and historic windows. Compatibility matrices are not associatedwith precise trajectories as are the cost function values. Compatibilitymatrices represent consistency or compatibility measures for anassignment of a sensor report to one trajectory in one frame with otherpossible assignments of other sensor reports in other frames. Suchmeasures apply to the construction of the set of all trajectoriesincorporating these linked sensor report assignments. In other words,individual components of any compatibility matrix are not in one-to-onecorrespondence with the set of possible trajectories as are thecomponents of the cost function. The system of varied-ordercompatibility matrices provides an extensive, flexible and varied datastructure for incorporating prior knowledge into the problem.

Third, Objective Function: The objective function which is to bemaximized is the average local consistency

$A = {{\sum\limits_{i,\Omega,\alpha}{{S\left( {i,\Omega,\alpha} \right)} \cdot {P\left( {i,\Omega,\alpha} \right)}}} = {{\sum\limits_{i,\Omega,\alpha}{\sum\limits_{j,\Gamma,\beta}{\ldots{\sum\limits_{Z}{{W_{Z} \cdot {R_{Z}\left( {i,\Omega,{\alpha{{j,\Gamma,\beta}}\ldots}} \right)} \cdot {P\left( {i,\Omega,\alpha} \right)} \cdot {P\left( {j,\Gamma,\beta} \right)}}{\ldots A}}}}}} = {\sum\limits_{i,\Omega,\alpha}{P_{i\Omega\alpha} \cdot S_{i\Omega\alpha}}}}}$

which measures the compliance of the real valued non-binary weightedlabeling assignments P(i,Ω,α) with the consistency conditionsrepresented by the system of compatibility matrices. The resultantS(i,Ω,α) support function is not a cost function analogous to C_(ijk) .. . but averages over the system of compatibility matrices with theweighted labeling assignments, P(i,Ω,α). Actually A is a nonlinearfunction of the weighted labeling assignments. Note that there is nosystem of constraints to complicate the relaxation labeling processes.

Fourth, once the system of compatibility matrices is formulated toinclude prior knowledge, simple updating formulae of the relaxationlabeling algorithm may be utilized without any need to recast theproblem to facilitate relaxation the relaxation process. Therefore, therelaxation labeling updating equations are invariant. This is in sharpcontrast to the complicating necessity of the Lagrangian relaxationtechniques to transform the problem into continuous sub-problems toeffect optimization. Also, simple conditions on formulating the systemof compatibility matrices result in very simple updating equations thatinsure non-negative changes to the objective function A and alsoguarantee the invariance of probability axioms,0≦P ^(K+1)(i,Ω,α)=P ^(K)(i,Ω,α)+dP ^(K)(i,Ω,α)≦1.0

It is the values of the weighted labeling assignments that are relaxedby the simpler relaxation labeling algorithm. There is a significantdifference between relaxing a set of values for P(i,Ω,α) with a simpleset of updating equations and relaxing a set of constraint equationsnecessitating a recasting of the problem formulation.

Fifth, Outputs: The real valued weighted labeling assignments P(i,Ω,α)which may or may not achieve a binary set of values necessary to defineprecise trajectories, i.e.

P(i,Ω,α)ε{0,1} and

${\sum\limits_{\alpha}{P\left( {i,\Omega,\alpha} \right)}} = {1.0.}$In general, under the simple relaxation labeling technique, the weightedlabeling assignments should approach close enough to this ideal to beuseful without change or to be easily adjusted to make them unambiguous.

Fifth, Coupling with ATR systems: Note that the formulation of thesystem of compatibility matrices provides facility for dynamicallycoupling the relaxation labeling processing and the ATR processing. Thisis why the ATR in ATRACK refers to Automatic Target Recognition and theTRACK refers to tracking. A cost function may incorporate inputs from anATR system but all ATR information pertaining to a track is factoredinto a single value of cost. In ATRACK, ATR information is factored intocompatibility matrices which link a small set of frames and notassociated with a trajectory which links all frames in a moving window.Most basic compatibility matrices are second order linking two reportsin two frames not necessarily sequential or adjacent frames.

ATRACK deals with fuzzy track objects because in every frame in theactive window each report may be associated with different weights withall trajectory indices and in each frame each trajectory object (index)may be associated with all reports and/or even labeled with the missingreport symbol. The weighted labeling assignments P(i,Ω,α) give therelative probability or weight of associating report a in frame Ω withtrack object i. The compatibility matrices give a consistency measurebetween two or more assignments and provide the means for incorporatingprior knowledge into the relaxation labeling approach. Differentcompatibility matrices can be conveniently used to incorporate differentaspects of prior knowledge including information provided by an ATRsystem. Factoring the compatibility matrices into linking matrices andATR dependent conditional probability distributions further facilitatesthe inclusion of ATR information. ATR information may also beincorporated into the linking matrices which are not associated withtrack objects. In ATRACK the conditional probability distribution is areport linking entity that factors into P(i,Ω,α) weighted averages ofthe identification vector P(c|Ω,α) and is therefore useful feedback datafor the ATR systems updating of subsequent identification vectors forATRCK to use in updating the conditional probability distribution.

In Poore's approach everything is factored into one cost value for eachpossible non-ambiguous trajectory. Further more it is not evident how toestablish a two way dynamic coupling between the Lagrangian relaxationprocess and the ATR updating process. If trajectory cost values areupdated with new ATR information the constraint equations are notchanged; therefore, the solution state must be checked for feasibilityand adjust to the new values of the cost function. The interestingquestion is how does the Lagrangian relaxation process provideintermediate solution states as useful feedback to the ATR processingmodule in order to affect a two way dynamic coupling process. The answerto this question is not discussed by Poore nor is evident how toduplicate ATRACK's dynamic coupling potential in which the optimizationof both ATR and track descriptors mutually assist one anotheriteratively

Superficial similarities exist between generally-known RLA and theinventive LRA due to nature of basic data formats for sensor outputs.First, sensor reports are the natural outputs of sensors. In ATRACTsensor reports also include data from an automatic target recognition(ATR) system which also processes the initial sensor reports.

Second, sensors generally make scans resulting in spatially andtemporally related reports. This naturally related group of sensorreports is called a frame. Even if sensor reports were gathered in acompletely random and unsystematic manner, it would be natural toarrange observed sensor reports into this quantized format.

Third, adjacent frames are also grouped and processed as active setswhich naturally defines a temporally moving window of active frames asolder frames are fixed and new frames acquired as active input. Onecould possibly work with an expanding window of frames but this wouldsoon become untenable with a continuous inflow of sensor reports. Thenumber of frames per active window is not generally a hard requirementbut usually adopted as a convenient simplification.

Fourth, therefore, these natural data groupings would be hard to avoidin any tracking problem definition which attempts to temporally linksensor reports initiated by the same sensed object. Utilization of thesedata structures does not limit the countless ways to incorporate priorknowledge, establish objectives influencing the track definitions,methods for optimizing selected objective function, and the type ofsolution achieved by the adopted processing technique.

As a further example, it will be appreciated that the RLA loopadvantageously allows a simple approach to optimizing the objectivefunction by optimizing the gradient and average of the gradient of theobjective function. However, other optimization approaches may beemployed.

1. A method of determining trajectories of a plurality of periodicallysensed targets, comprising: periodically accepting a plurality of targetsensor reports (α, β, χ, . . . ); grouping temporally related subsets ofthe plurality of target sensor reports into a plurality a sequentialframes (Ω, Γ, . . . ); establishing a system of track objects (i, j, k,. . . ) associated with the plurality of target sensor reports in aprevious frame; linking each target sensor report (α) in the currentframe to each of one or more track objects {i, j, . . . } in a previousframe (Ω) by a weighted labeling assignment (P(i, Ω, α)); defining aplurality (M) of compatibility matrices (R_(Z)(i, Ω, α|j, Γ, β| . . . ))wherein Z=1,2, . . . ,M) each operatively configured to measure selectedrelative compatibilities and consistencies of a set of sensor reportslabels identified by {(Ω, α), (Γ, β), . . . } being assigned to the setof track objects identified by {i, j, . . . }; generating a system ofsupport functions S_(z)(i, Ω, α) from said system of varied-orderedcompatibility matrices (R_(z)) by computing weighted averages of eachcompatibility matrix (R_(z)) with said weighted labeling assignments,$\left( {{S_{z}\left( {i,\Omega,\alpha} \right)} = {\sum\limits_{j,\Gamma,\beta}{\ldots\mspace{14mu}{{R_{z}\left( {i,\Omega,{\alpha{{j,\Gamma,\beta}}\ldots}} \right)} \cdot {P\left( {j,\Gamma,\beta} \right)}}\ldots}}} \right)$generating a resultant support function [S(i,Ω,α)] by combining saidsystem of support functions S_(z) (i,Ω,α); iteratively modifying saidweighted labeling assignments (P(i,Ω,α)) with a relaxation labelingtechnique that maximizes an average local consistency$\left( {A = {\sum\limits_{i,\Omega,\alpha}{{P\left( {i,\Omega,\alpha} \right)} \cdot {S\left( {i,\Omega,\alpha} \right)}}}} \right)$which measures compliance of said weighted labeling assignments withsaid consistency conditions represented by said system of compatibilitymatrices; generating trajectory descriptors of periodically sensedtargets using said weighted labeling assignments.
 2. The method of claim1, wherein defining the plurality (M) of compatibility matrices(R_(z)(i, Ω, α|j, Γ, β| . . . ) where Z=1,2 . . . ,M) iterativelymodifying said weighted labeling assignments P(i, Ω, α) with therelaxation labeling technique that maximizes an average localconsistency further comprises: formulating said system of varied-orderedcompatibility matrices having symmetric structure and complying withconditions−1.0≦L ^(Z) _(MIN) ≦R _(Z)(i,Ω,α|j,Γ,β . . . )≦L ^(Z) _(MAX)≦1.0 for allZ,i,Ω,α,j,Γ,β, . . . ,  where L^(Z) _(MAX)−L^(Z) _(MIN)≦1.0 and Z=1.2, .. . ,M; initializing said weighted labeling assignments to satisfyprobability axioms; i.e., 0≦P(i,Ω,α)≦1.0 and${{\sum\limits_{\alpha}{P\left( {i,\Omega,\alpha} \right)}} = 1.0};$generating said resultant support function by combining said system ofsupport functions by weighted average; i.e.,${{S^{K}\left( {i,\Omega,\alpha} \right)} \equiv {\sum\limits_{Z}{W_{Z} \cdot {S_{Z}^{K}\left( {i,\Omega,\alpha} \right)}}}},$where W_(Z) satisfies said probability axioms and K=iteration index;generating averages of said resultant support function,${{{\hat{S}}^{K}\left( {i,\Omega} \right)} \equiv {\sum\limits_{\alpha}{{S^{K}\left( {i,\Omega,\alpha} \right)} \cdot {P^{K}\left( {i,\Omega,\alpha} \right)}}}};$generating modifications to said weighted labeling assignments with saidrelaxation labeling technique defined bydP^(K)(i,Ω,α)=P^(K)(i,Ω,α)·(G^(K)(i,Ω,α)−Ĝ^(K)(i,Ω))${{G^{K}\left( {i,\Omega,\alpha} \right)} = {{\frac{\partial A}{\partial{P^{K}\left( {i,\Omega,\alpha} \right)}} \cdot \frac{1}{\sum\limits_{Z}{W_{Z} \cdot O_{Z}}}} = {\sum\limits_{Z}{V_{Z} \cdot {S_{Z}^{K}\left( {i,\Omega,\alpha} \right)}}}}},{{{\hat{G}}^{K}\left( {i,\Omega} \right)} = {\sum\limits_{\alpha}{{P^{K}\left( {i,\Omega,\alpha} \right)} \cdot {G^{K}\left( {i,\Omega,\alpha} \right)}}}},{where}$${V_{Z} = {{\frac{W_{Z} \cdot O_{z}}{\sum\limits_{Z}{W_{Z} \cdot O_{Z}}}{and}\mspace{14mu} O_{z}} = {{Order}\mspace{14mu}{of}\mspace{14mu} R_{Z}}}};$updating said weighted labeling assignments with said modificationsdP^(K)(i,Ω,α) which automatically preserves said initialing probabilityaxioms; i.e., 0≦P^(K+1)(i,Ω,α)=P^(K)(i,Ω,α)+dP^(K)(i,Ω,α)≦1.0 and${{\sum\limits_{\alpha}{P^{K + 1}\left( {i,\Omega,\alpha} \right)}} = 1.0},$and automatically insures non-negative changes to said average localconsistency measure; i.e., $\begin{matrix}{A = {\sum\limits_{i,\Omega,\alpha}{{S^{K}\left( {i,\Omega,\alpha} \right)} \cdot {P^{K}\left( {i,\Omega,\alpha} \right)}}}} \\{{= {\sum\limits_{i,\Omega,\alpha}{\sum\limits_{j,\Gamma,\beta}{\ldots{\sum\limits_{Z}{{W_{Z} \cdot {R_{Z}\left( {i,\Omega,{\alpha{{j,\Gamma,\beta}}\ldots}} \right)} \cdot {P^{K}\left( {i,\Omega,\alpha} \right)} \cdot {P^{K}\left( {j,\Gamma,\beta} \right)}}\ldots}}}}}}\mspace{14mu},}\end{matrix}$${{dA}{\sum\limits_{i,\Omega,\alpha}{\frac{\partial A}{\partial{P^{K}\left( {i,\Omega,\alpha} \right)}} \cdot {{dP}^{K}\left( {i,\Omega,\alpha} \right)}}}} \geq {0.0.}$3. The method of claim 1, wherein linking each report in the frame toone or more tracks in a previous frame further comprises performingautomatic target recognition analysis to generate a probability ofassignment for each report.
 4. The method of claim 1, wherein linkingeach report in the frame to one or more tracks in a previous framefurther comprises performing automatic target recognition analysis togenerate a probability of assignment for each report.
 5. The method ofclaim 4, wherein performing automatic target recognition analysiscomprises performing target signature correlation.
 6. The method ofclaim 4, wherein performing automatic target recognition analysiscomprises performing predictive matching.
 7. The method of claim 4,wherein performing automatic target recognition analysis comprisesperforming mobility analysis with reference to a characteristicsdatabase.
 8. The method of claim 4, wherein performing automatic targetrecognition analysis comprises performing co-moving target analysis withreference to a characteristics database.
 9. The method of claim 1,wherein linking each report in the frame to one or more tracks in aprevious frame further comprises estimating masking areas with referenceto a terrain database and computing a compatibility matrix for theestimated masking areas.
 10. The method of claim 1, wherein convertingthe linking matrix into the plurality of compatibility matrices furthercomprises computing a negative compatibility matrix for encouragingcompetition between assignment vectors to a given report.
 11. The methodof claim 1, wherein converting the linking matrix into the plurality ofcompatibility matrices further comprises coupling a targetidentification vector with the assignment vector.
 12. The method ofclaim 1, wherein defining said system of varied-ordered compatibilitymatrices to incorporate knowledge of consistency conditions influencinglabeling of said track objects further comprises: accepting inputs froman automatic target recognition (ATR) system; and defining a system ofvaried-ordered compatibility matrices incorporating said ATR systeminputs.
 13. The method of claim 12, wherein forming said system ofcompatibility matrices combining inputs from said automatic targetrecognition (ATR) system further comprises: establishing a system ofvaried-ordered lining matrices independent of said track objects;combining said varied-ordered linking matrices with said inputs fromsaid ATR system forming a corresponding subsystem of said varied-orderedATR dependent compatibility matrices.
 14. The method of claim 13,wherein establishing a system of varied-ordered linking matrices furthercomprises: defining a system of varied-ordered linking matricesincorporating selected aspects of prior knowledge with respect torelative consistency of linking said sensor reports for target classesidentifiable by class index c; and combining said system of varied-orderlinking matrices into a system of resultant linking matrices, whereL_(N)(Ω,α|Λ,β| . . . |c) identifies said resultant linking matrix fororder N and target class c.
 15. The method of claim 13, whereingenerating said system of varied-ordered compatibility matrices bycombining said system of linking matrices with said automatic targetrecognition (ATR) system inputs further comprises: accepting said ATRsystem inputs; formulating with said ATR system inputs a conditionalprobability distribution designated by P(c|j); convert said resultantlinking matrix of varied-order N into said AIR dependent compatibilitymatrix of varied-order N utilizing said conditional probabilitydistribution,${R_{N}\left( {j,\Omega,{\alpha{{j,\Lambda,\beta}}\ldots}} \right)} \equiv {\sum\limits_{c}{{L_{N}\left( {\Omega,{\alpha{{\Lambda,\beta}}\ldots{\left. c \right) \cdot {P\left( c \right.}}j}} \right)}.}}$16. The method of claim 15, wherein accepting said automatic targetrecognition (ATR) system inputs and formulating said conditionalprobability distribution further comprises: accepting as part of saidATR system inputs an ATR system included identification vectordesignated by P(c|Ω,α); and defining averaged conditional probabilitydistributions identifiable with <P(c|j)> by taking averages of saididentification vector P(c|Ω,α) with said weighting labeling assignmentsP(i,Ω,α).
 17. The method of claim 13, wherein generating said system ofvaried-ordered compatibility matrices by combining said system ofvaried-ordered linking matrices with said automatic target recognition(ATR) system inputs further comprises: periodically updating saidaveraged conditional probability distributions combining iterativeenhancements to said weighting labeling assignments resulting from saidrelaxation labeling processing with periodic ATR generated enhancementsto said ATR system inputs; and applying said enhanced averagedconditional probability distributions to updating said system of AIRdependent varied-ordered compatibility matrices.
 18. The method of claim13, wherein generating said system of varied-ordered compatibilitymatrices by combining said system of varied-ordered linking matriceswith said automatic target recognition (ATR) system inputs furthercomprises: storing said varied-ordered linking matrices instead of saiddependent compatibility matrices resulting in reduced storagerequirements; and computing said support functions utilizingcorresponding said stored linking matrices.
 19. The method of claim 12,wherein accepting said inputs from said automatic target recognition(ATR) system and defining said system of varied-ordered compatibilitymatrices incorporating said ATR inputs further comprises performingautomatic target recognition analysis to generate relative measures forlinking two or more said sensor reports.
 20. The method of claim 19,wherein performing automatic target recognition analysis furthercomprises performing target signature correlation.
 21. The method ofclaim 19, wherein performing automatic target recognition analysisfurther comprises performing predictive matching.
 22. The method ofclaim 19, wherein performing automatic target recognition analysisfurther comprises performing mobility analysis with reference to acharacteristics database.
 23. The method of claim 19, wherein performingautomatic target recognition analysis further comprises performingco-moving target analysis with reference to a characteristics database.24. The method of claim 1, wherein generating trajectory descriptors ofperiodically sensed targets using said weighted labeling assignmentsmodified with said relaxation labeling technique and outputting saiddetermined trajectories descriptors further comprises: accepting inputsfrom an automatic target recognition (ATR) system to enhance saidtrajectory descriptors; and outputting said trajectories descriptorsincluding said updated averaged conditional probability distributions tosaid ATR system for enhancing subsequent ATR generated inputs, thuscoupling said ATR system processing and said relation labelingprocessing and reciprocally enhancing results from both processingsystems.
 25. The method of claim 1, wherein defining higher ordercompatibility matrices further comprises estimating marking areas withreference to a characteristics database to account for consistency ofsaid sensor and missing sensor reports and computing a compatibilitymatrix for said estimated masking areas.
 26. The method of claim 1,wherein defining higher order compatibility matrices further comprisesconstructing compatibility matrices combining other said compatibilitymatrices of lower-order as defining components.
 27. The method of claim1, wherein defining said system of varied-order compatibility matricesfurther comprises adjusting maximum separation between frames linked bysaid system of multi-ordered compatibility matrices.
 28. The method ofclaim 1, wherein iteratively modifying said weighted labelingassignments with said relaxation labeling technique further comprises:restricting said iterative modification of said weighted labelingassignments to a set of active components of said weighted labelingassignments associated with a temporally changing set of active framesdesignated as an active window; restricting said iterative modificationof said weighted labeling assignments in inactive frames outside saidactive window referred to as an inactive window.
 29. The method of claim1, wherein defining said system of varied-order compatibility matricesfurther comprises redefining a set of said compatibility matricesconnecting frames in both said active and inactive windows as a set ofreduced compatibility matrices having reduced order equaling number ofsaid active frames connected.
 30. The method of claim 1, whereingenerating a system of support functions from said system ofvaried-ordered compatibility matrices further comprises: maintaininginactive components of said weighted labeling assignments associatedwith said inactive window for use in generating support functions; andcomputing support functions from said set of reduced multi-ordercompatibility matrices for defining said resultant support function useby said relaxation labeling processing in said active window.
 31. Themethod of claim 1, wherein establishing system of track objects dryercomprises methods for initiating said track objects consistent with saidsensor reports.
 32. The method of claim 1, wherein establishing saidsystem of tracking objects further comprises methods for removal of saidtrack objects no longer consistent with said sensor reports.
 33. Themethod of claim 1, wherein defining said system of compatibilitymatrices further comprises initializing new components of saidcompatibility matrices associated with said initializing of trackobjects.
 34. The method of claim 1, wherein defining said system ofcompatibility matrices further comprises initializing values for saidcompatibility matrix components associated with said missing sensorreport labels.
 35. The method of claim 1, wherein defining weightedlabeling assignments further comprises initializing weighted labelingassignments associated with said initializing of new components of saidcompatibility matrices.
 36. The method of claim 1, wherein defining saidsystem of compatibility matrices further comprises dynamically alteringa number of varied-ordered compatibility matrices defining saidresultant support function.
 37. The method of claim 1, wherein definingsaid system of varied-ordered compatibility matrices further comprisesdynamically altering definitions and orders of said compatibilitymatrices defining said resultant support function.
 38. The method ofclaim 1, wherein fixing values of said weighted labeling assignments todetermine said trajectories further comprises adjusting said fixedvalues to selected track-defining levels of precision ranging fromallowing strictly ambiguous values for fuzzy trajectory definition to aset of binary values for precise unambiguous trajectory definitiondepending on requirements of said external application modules utilizingsaid trajectory definitions as input.
 39. An apparatus, comprising: amemory; a plurality of target sensor reports (α, β, χ, . . . ) residentin the memory; and a program resident in the memory, the programconfigured (a) to group temporally related subsets of the plurality oftarget sensor reports into a plurality a sequential frames (Ω, Γ, . . .); (b) to establish a system of track objects (i, j, k, . . . )associated with the plurality of target sensor reports in a previousframe; (c) to link each target sensor report (c) in the current frame toeach of one or more track objects (i) in a previous frame (Ω) by aweighted labeling assignment P(i, Ω, α)): (d) to define a plurality (M)of compatibility matrices (R_(Z)(i, Ω, α|j, Γ, β| . . . ) where Z=1,2, .. . ,M) operatively configured to measure selected relativecompatibilities and consistencies of a set of sensor reports labelsidentified by {(Ω, α), (Γ, β), . . . } being assigned to a set of trackobjects identified by {i, j, . . . }; (e) to generate a system ofsupport functions S_(Z)(i, Ω, α) from said system of varied-orderedcompatibility matrices (R_(Z)) by computing weighted averages of eachcompatibility matrix (R_(Z)) with said weighted labeling assignments.$\left( {{S_{z}\left( {i,\Omega,\alpha} \right)} = {\sum\limits_{j,\Gamma,\beta}{\ldots\mspace{14mu}{{R_{z}\left( {i,\Omega,{\alpha{{j,\Gamma,\beta}}\ldots}} \right)} \cdot {P\left( {j,\Gamma,\beta} \right)}}\ldots}}} \right)$generating a resultant support function (S(i,Ω,α)) by combining saidsystem of support functions S_(Z)(i,Ω,α); (f) to iteratively modify saidweighted labeling assignments (P(i, Ω, α)) with a relaxation labelingtechnique that maximizes an average local consistency$\left( {A = {\sum\limits_{i,\Omega,\alpha}{{P\left( {i,\Omega,\alpha} \right)} \cdot {S\left( {i,\Omega,\alpha} \right)}}}} \right)$which measures compliance of said weighted labeling assignments withsaid consistency conditions represented by said system of compatibilitymatrices; and (g) to generate trajectory descriptors of periodicallysensed targets using said weighted labeling assignments.
 40. A programproduct, comprising: (a) a program configured to (a) to group temporallyrelated subsets of a plurality of target sensor reports into a pluralitya sequential frames (Ω, Γ, . . . ); (b) to establish a system of trackobjects (i, j, k, . . . ) associated with the plurality of target sensorreports in a previous frame; (c) to link each target sensor report (α)in the current frame to each of one or more track objects (i) in aprevious frame (Ω) by a weighted labeling assignment (P(i, Ω, α)); (d)to define a plurality (M) of compatibility matrices (R_(Z)(i, Ω, α|j, Γ,β, . . . ) where Z=1,2, . . . ,M) operatively configured to measureselected relative compatibilities and consistencies of a set of sensorreports labels identified by {(Ω, α), (Γ, β), . . . } being assigned toa set of track objects identified by {i, j, . . . }; (e) to generate asystem of support functions S_(Z)(i, Ω, α) from said system ofvaried-ordered compatibility matrices (R_(Z)) by computing weightedaverages of each compatibility matrix (R_(Z)) with said weighted labelassignments,$\left( {{S_{z}\left( {i,\Omega,\alpha} \right)} = {\sum\limits_{j,\Gamma,\beta}{\ldots\mspace{14mu}{{R_{z}\left( {i,\Omega,{\alpha{{j,\Gamma,\beta}}\ldots}} \right)} \cdot {P\left( {j,\Gamma,\beta} \right)}}\ldots}}} \right)$generating a resultant support function (S(i,Ω,α)) by combining saidsystem of support functions S_(Z)(i,Ω,α); (f) to iteratively modify saidweighted labeling assignments (P(i,Ω,α)) with a relaxation labelingtechnique that maximizes an average local consistency$\left( {A = {\sum\limits_{i,\Omega,\alpha}{{P\left( {i,\Omega,\alpha} \right)} \cdot {S\left( {i,\Omega,\alpha} \right)}}}} \right)$which measures compliance of said weighted labeling assignments withsaid consistency conditions represented by said system of compatibilitymatrices; and (g) to generate trajectory descriptors of periodicallysensed targets using said weighted labeling assignments; and (b) asignal bearing media bearing the program.
 41. The program product ofclaim 40, wherein the signal bearing media comprises at least one of arecordable media and a transmission-type media.